sumt_2.pas

机械优化设计中的惩罚函数法

unit sumt_2;  //SUMT_算法主程序

interface

uses
  Windows, Messages, SysUtils, Classes, Graphics, Controls, Forms, Dialogs,
  StdCtrls, ComCtrls,sumt_1, sumt_0;

type
  TForm2 = class(TForm)
    GroupBox1: TGroupBox;
    cp: TButton;
    qx: TButton;
    sd: TSaveDialog;
    jgxs: TMemo;
    procedure FormCreate(Sender: TObject);
    procedure cpClick(Sender: TObject);
    procedure qxClick(Sender: TObject);
  private
    { Private declarations }
  public
    { Public declarations }
  end;

var
  Form2: TForm2;
  procedure sjxs_1;
  procedure sjxs_2;
  Procedure sump;
  procedure xran;
  procedure powell;
  procedure dfpbfg;
  procedure rvalul;
  procedure linmin;
  procedure functi;
  procedure grapen;
  procedure differ;
  function rmd() : real;
  function abc1(n:integer;a,b:sumt_1.arr1):real;
  function abc2(ii,n:integer;H:sumt_1.arr2;a:sumt_1.arr1):real;

implementation

uses sumt_fgh;

{$R *.DFM}

procedure TForm2.FormCreate(Sender: TObject);  //第二幅窗体启动
begin
   Form2.Left :=50;
   Form2.Top  :=50;
end;

procedure TForm2.cpClick(Sender: TObject);    //存盘
begin
  if sd.Execute then
  begin
    if fileexists(sd.filename) then
    begin
      if application.MessageBox('您确认要覆盖此文件吗？','警告',MB_YESNO)=idyes then
        form1.WriteDataToFile(sd.filename);
    end
    else
    begin
      form1.WriteDataToFile(sd.filename);
    end;
    if application.MessageBox('要返回吗？','提示',MB_YESNO)=idyes then Form2.Close;
  end;
end;

procedure TForm2.qxClick(Sender: TObject);  //取消或返回
begin
//  form2.Destroy;
//  Application.CreateForm(TForm2, Form2);
  Form2.Close;
end;

//=================优化算法程序=============================
procedure sump;   //算法主辅程序
var i,k,temp:integer;xx,yy,pen0:real;
    x0:arr1;
    label Here;
begin
  with form1.sumt do
  begin
    ffx;ggx;
    ffxx0:=fx;
    temp:=0;
    for k:=1 to kg do if gx[k]>=0.0 then begin xran; break; end;
    if (R<=0.0) then Rvalul; functi;
    for i:=1 to n do x0[i]:=x[i];  //pen0:=pen;
    with form2.jgxs.lines do
    begin
      add('-------------------------------------------------------------------------------');
      add('                                                              ');
      add('二、计算过程__数据');
      add('===============================================================================');
    repeat
      for i:=1 to n do x0[i]:=x[i];
      pen0:=pen;

      Add('   IRC = '+inttostr(IRC)+#9+
          '   R = '+formatfloat('0.000000E+00',R)+#9+
          '   PEN = '+floattostr(PEN));
      Add('  -----------------------------------------------------------------------------');

      if kcheck=0 then powell
                  else dfpbfg;
      ITE:=ITE+KTE;
      r:=r*cr;
      IRC:=IRC+1;
      for i:=1 to n do
      begin
        xx:=abs(x0[i]-x[i]); yy:=abs(x[i])*eps+eps;
        if xx>yy then  begin temp:=1; break ;end
                 else  temp:=0;
      end;
    until (abs(pen-pen0) <= eps*abs(pen)+eps) or (temp=0);
  end;
  end;
end;

function rmd() : real;  //产生随机数程序段
var
   rm35,rm36,rm37 :real;
begin
  with form1.sumt do
  begin
    rm35:=exp(35.0*ln(2.0));
    rm36:=2.0*rm35;
    rm37:=2.0*rm36;
    rm:=5.0*rm;
    if rm>=rm37 then rm:=rm-rm37;
    if rm>=rm36 then rm:=rm-rm36;
    if rm>=rm35 then rm:=rm-rm35;
    rmd:=rm/rm35;
  end;
end;

function ppp(aa,bb:real) : real;
var cc:real;
begin
 cc:=bb*ln(aa);
 ppp:=exp(cc);
end;

procedure xran;   //随机产生可行的初始点;
var ii,temp1,kk:integer;qq:real;
begin
  with form1.sumt do
  begin
    with form2.jgxs.lines do
    begin
      add('###############################################################################');
      add('  初始点不可行！');
      Add('  -----------------------------------------------------------------------------');
      add('  不可行的初始点：');
      Add('  -----------------------------------------------------------------------------');
      for ii := 1 to n do Add(#9+#9+#9+'X['+inttostr(ii)+']='+floattostr(x[ii]));
      Add('  -----------------------------------------------------------------------------');
      Add('   初始点处的不等约束函数值 G(X*):');
      for ii := 1 to kg do Add(#9+#9+#9+#9+'GX['+inttostr(ii)+']= '+formatfloat('0.000000E+00',GX[ii]));
    end;

    repeat
      temp1:=0;
      ggx;
      for ii:=1 to kg do
      begin
        if (gx[ii]>=0) then
        begin  for kk:=1 to n do
               begin  qq:=rmd();  x[kk]:=bl[kk]+qq*(bu[kk]-bl[kk]); end;
               temp1:=1;
               break;
        end;
      end;
    until temp1=0;

    with form2.jgxs.lines do
    begin
      Add('  -----------------------------------------------------------------------------');
      add('  找到可行的初始点！');
      Add('  -----------------------------------------------------------------------------');
      add('  可行的初始点：');
      Add('  -----------------------------------------------------------------------------');
      for ii := 1 to n do Add(#9+#9+#9+'X['+inttostr(ii)+']='+floattostr(x[ii]));
      Add('  -----------------------------------------------------------------------------');
      Add('   初始点处的不等约束函数值 G(X*):');
      for ii := 1 to kg do Add(#9+#9+#9+#9+'GX['+inttostr(ii)+']= '+formatfloat('0.000000E+00',GX[ii]));
      add('###############################################################################');
      add(' ');
    end;

  end;
end;

procedure rvalul;    //计算初始惩罚因子程序段
var i:integer;pin:real;
begin
  with form1.sumt do begin
    ffx;   ggx;
    pin:=0.0;
    for i:=1 to kg do pin:=pin+1.0/abs(gx[i]);
    R:=abs(fx/pin);
  end;
end;

procedure functi;   //计算惩罚函数值程序段
var i:integer;  pin,pex:real;
begin
  with form1.sumt do
  begin
    NPE:=NPE+1;
    ffx;  ggx;  hhx;
    pin:=0.0;  pex:=0.0;
    if kg>0 then for i:=1 to kg do pin:=pin+1.0/gx[i];
    if kh>0 then for i:=1 to kh do pex:=pex+hx[i]*hx[i];
    pen:=fx-R*pin+pex/sqrt(R);
  end;
end;

procedure grapen;  //计算惩罚函数对X的偏导数
var i,k:integer;  din,dex,dfu:real;
begin
  with form1.sumt do
  begin
    differ;
    for i:=1 to n do
    begin
      din:=0.0; dex:=0.0; dfu:=dfdx[i];
      for k:=1 to kg do din:=din-dgdx[i,k]/(gx[k]*gx[k]);
      for k:=1 to kh do dex:=dex+2.0*dhdx[i,k]*hx[k];
      dpdx[i]:=dfu+r*din+dex/sqrt(R);
    end;
  end;
end;

procedure differ;  //用中心差商法求偏导数
var i,k:integer;dj,xi:real;
begin
  with form1.sumt do
  begin
    ffx;  ggx;  hhx;
    for i:=1 to n do
    begin
      dj:=1e-6*abs(x[i]);
      if dj>abs(fx) then dj:=abs(fx);
      for k:=1 to kg do if dj>abs(gx[k]) then dj:=abs(gx[k]);
      for k:=1 to kh do if dj>abs(hx[k]) then dj:=abs(hx[k]);
      if dj<1e-5 then dj:=1e-5;
      xi:=x[i];x[i]:=xi+dj;ffx;ggx;hhx;
      dfdx[i]:=fx;
      for k:=1 to kg do dgdx[i,k]:=gx[k];
      for k:=1 to kh do dhdx[i,k]:=hx[k];
      x[i]:=xi-dj;ffx;ggx;hhx;
      dfdx[i]:=(dfdx[i]-fx)/(2.0*dj);
      for k:=1 to kg do dgdx[i,k]:=(dgdx[i,k]-gx[k])/(2.0*dj);
      for k:=1 to kh do dhdx[i,k]:=(dhdx[i,k]-hx[k])/(2.0*dj);
      x[i]:=xi;
    end;
  end;
end;

procedure linmin;   //一维搜索程序段
var i,k,ch,temp1:integer;
    T1,T2,T3,T4,TH,FF1,FF2,FF3,FF4,C1,C2:real;
    xx:arr1;
    label Here1;
begin
  with form1.sumt do
  begin
    ILI:=ILI+1;
    temp1:=0;
    T1:=0.0; TH:=T0; T2:=T0; FF1:=pen;
    for i:=1 to n do xx[i]:=x[i];
    repeat
      for i:=1 to n do x[i]:=xx[i]+T2*S[i];
      ggx;
      for k:=1 to kg do
        if gx[k]>0.0 then  begin  temp1:=1;T2:=0.5*T2;break; end
                     else temp1:=0;
    until temp1=0;
    functi;  FF2:=pen;
    if FF2>FF1 then   //换向
    begin
        TH:=-TH;
        T3:=T1; FF3:=FF1;    T1:=T2; FF1:=FF2;   T2:=T3; FF2:=FF3;
    end;
    T3:=T2+TH;
    repeat
      for i:=1 to n do x[i]:=xx[i]+T3*S[i];
      ggx;
      for k:=1 to kg do
        if gx[k]>0.0 then
         begin temp1:=1;T3:=T3-0.5*(T3-T2);
           if (abs(T3-T2)<=1e-7) then
            begin
              for i:=1 to n do x[i]:=xx[i]+T2*S[i];pen:=FF2; goto here1;
            end;
            break;
         end else temp1:=0;
    until temp1=0;
    functi;FF3:=pen;
    while ff2>ff3 do begin  //加速、跟进、外推
      TH:=TH+TH;   T1:=T2; FF1:=FF2;  T2:=T3; FF2:=FF3;   T3:=T2+TH;
    repeat
      for i:=1 to n do x[i]:=xx[i]+T3*S[i];
      ggx;
      for k:=1 to kg do
        if gx[k]>0.0 then
        begin
          temp1:=1;T3:=T3-0.5*(T3-T2);
          if (abs(T3-T2)<=1e-7) then
          begin
            for i:=1 to n do x[i]:=xx[i]+T2*S[i];
            pen:=FF2;
            goto here1;
          end;
          break;
        end
        else temp1:=0;
      until temp1=0;
      functi;   FF3:=pen;
    end;
    repeat       //以下为二次插值法程序
      C1:=(FF3-FF1)/(T3-T1);
      if (abs(T2-T3)<=1e-7) then
      begin
         for i:=1 to n do x[i]:=xx[i]+T2*S[i];
         pen:=FF2;
         goto here1;
      end;
      C2:=((FF2-FF1)/(T2-T1)-C1)/(T2-T3);
      if (abs(C2)<=1e-7) then
      begin
        for i:=1 to n do x[i]:=xx[i]+T2*S[i];
        pen:=FF2;
        goto here1;
      end;
      T4:=0.5*(T1+T3-C1/C2);
      if ((T4-T1)*(T3-T4)<=0.0) then
      for i:=1 to n do
      begin