📄 sumt_2.pas
字号:
x[i]:=xx[i]+T2*S[i];
pen:=FF2;
goto here1;
end;
for i:=1 to n do x[i]:=xx[i]+T4*S[i];
functi; FF4:=pen;
begin
ch:=0;
if((T4-T2)*TH >0.0) and(FF2>FF4) then ch:=1;
if((T4-T2)*TH >0.0) and(FF2<=FF4)then ch:=2;
if((T4-T2)*TH<=0.0) and(FF2>FF4) then ch:=3;
if((T4-T2)*TH<=0.0) and(FF2<=FF4)then ch:=4;
case ch of
1: begin T1:=T2;FF1:=FF2;T2:=T4;FF2:=FF4; end;
2: begin T3:=T4;FF3:=FF4; end;
3: begin T3:=T2;FF3:=FF2;T2:=T4;FF2:=FF4; end;
4: begin T1:=T4;FF1:=FF4; end;
end;
end;
until (abs(FF2-FF4)<(eps*abs(FF2)+eps))or
(abs(T2-T1)<=1e-7) or (abs(T2-T3)<=1e-7);
if(FF4<FF2) then
for i:=1 to n do begin x[i]:=xx[i]+T4*S[i]; pen:=FF4; end
else
for i:=1 to n do begin x[i]:=xx[i]+T2*S[i]; pen:=FF2; end;
here1:
end;
end;
procedure powell; //POWELL法
var i,j,k,idf,n1,i1,temp1:integer;
af,bf,F,F1,F2,F3,dfm,df,sum:real;
x1,x2,x3,y:arr1;
aaa,bbb,cc1,ddd:real;
begin
with form1.sumt do
begin
KTE:=0;
temp1:=0;
for i:=1 to n do
begin
for k:=1 to n do direct[k,i]:=0.0;
direct[i,i]:=1.0;
end;
repeat
KTE:=KTE+1;
F :=pen; F1:=F; F2:=F;
for i:=1 to n do begin x1[i]:=x[i]; x2[i]:=x[i]; y[i]:=x[i];end;
idf:=1;
dfm:=0.0;
for i:=1 to n do
begin
for k:=1 to n do S[k]:=direct[k,i];
for j:=1 to n do x[j]:=x2[j];
linmin;
for j:=1 to n do x2[j]:=x[j];
ggx;
for j:=1 to kg do
if gx[j]>0.0 then begin temp1:=1;break; end
else temp1:=0;
F2:=pen;
df:=F-F2;
F:=F2;
if(df>dfm) then begin dfm:=df;idf:=i;end;
end;
for i:=1 to n do begin x3[i]:=2.0*x2[i]-x1[i]; x[i]:=x3[i];end;
ggx;
for k:=1 to kg do
if gx[k]>0.0 then begin temp1:=1;break;end
else temp1:=0;
if temp1=0 then
begin
functi;
F3:=pen;
af:=(F1-2.0*F2+F3)*(F1-F2-dfm)*(F1-F2-dfm);
bf:=0.5*dfm*(F1-F3)*(F1-F3);
if (F3<F1)and(af<bf) then
begin
n1:=n-1;
if(idf<>n) then
begin
for i:=idf to n1 do
begin
i1:=i+1;
for k:=1 to n do direct[k,i]:=direct[k,i1];
end;
end;
sum:=0.0;
for k:=1 to n do
begin
direct[k,n]:=x3[k]-x1[k];
sum:=sum+direct[k,n]*direct[k,n];
end;
sum:=1.0/sqrt(sum);
for k:=1 to n do direct[k,n]:=direct[k,n]*sum;
for k:=1 to n do s[k]:=direct[k,n];
for i:=1 to n do x[i]:=x2[i];
pen:=F2;
linmin;
for i:=1 to n do x2[i]:=x[i];
// F2:=pen;
end else
if(F3<F2) then
begin
for i:=1 to n do x[i]:=x3[i];
pen:=F3;
f:=f3;
end
else
begin
for i:=1 to n do x[i]:=x2[i];
pen:=F2;
f:=f2;
end;
end else
begin
for i:=1 to n do x[i]:=x2[i];
pen:=F2;
f:=f2;
end;
ffx; ggx; hhx;
for j:=1 to n do
if gx[j]>0.0 then begin temp1:=1;break;end else temp1:=0;
aaa:=abs(F-F1);
bbb:=abs(F)*eps+eps;
cc1:=0.0;
ddd:=0.0;
for i:=1 to kg do
begin
cc1:=sqrt(abs(x[i]*x[i]-y[i]*y[i]));
ddd:=sqrt(x[i]*x[i]);
end;
until (aaa<=bbb) or (cc1<=ddd);
end;
end;
procedure dfpbfg; //DFP法或BFGS法
var i,ii,k,ip:integer;F1, psi,a,sum:real;
label Here;
begin
with form1.sumt do
begin
KTE:=0;
here:
for i:=1 to n do
begin
for k:=1 to n do H[i,k]:=0.0;
H[i,i]:=1.0;
end;
grapen; F0:=pen; ip:=1;
repeat
KTE:=KTE+1;
if ip<>1 then
begin
for i:=1 to n do
begin dx[i]:=x[i]-dx[i];dg[i]:=dpdx[i]-dg[i];end;
if kcheck=1 then
begin
dxtdg:=abc1(n,dx,dg);
for ii:=1 to n do hdg[ii]:=abc2(ii,n,H,dg);
dgthdg:=abc1(n,dg,hdg);
for i:=1 to n do for k:=1 to n do
begin
H[i,k]:=H[i,k]+dx[i]*dx[k]/dxtdg-hdg[i]*hdg[k]/dgthdg;
if i<>k then H[k,i]:=H[i,k];
end;
end;
if kcheck=2 then
begin
for ii:=1 to n do hdg[ii]:=abc2(ii,n,H,dg);
dgthdg:=abc1(n,dg,hdg);
dxtdg:=abc1(n,dx,dg);
psi:=1.0+dgthdg/dxtdg;
for i:=1 to n do
for k:=1 to n do
begin
H[i,k]:=h[i,k]+((dx[i]*psi-hdg[i])*dx[k]-hdg[k]*dx[i])/dxtdg;
if i<>k then h[k,i]:=h[i,k];
end;
end;
end;
ip:=0;
F1:=pen;
for i:=1 to n do begin dx[i]:=x[i]; dg[i]:=dpdx[i];end;
for ii:=1 to n do s[ii]:=abc2(ii,n,H,dpdx);
sum:=0.0;
for i:=1 to n do sum:=sum+s[i]*s[i];
sum:=1.0/sqrt(sum);
a:=abc1(n,dpdx,s);
if a>0.0 then a:=-1.0 else a:=1.0;
for i:=1 to n do s[i]:=a*s[i]*sum;
linmin;
grapen;
if pen>F0 then goto here;
until (abs(pen-F1)<=eps*abs(pen)+eps);
end;
end;
function abc1(n:integer;a,b:arr1):real; //向量A 与向量B的数量积
var i:integer;
atb:real;
begin
atb:=0.0;
for i:=1 to n do atb:=atb+a[i]*b[i];
abc1:=atb;
end;
function abc2(ii,n:integer;H:arr2;A:arr1):real; //矩阵H与向量A的乘积
var j:integer;ha:real;//a:arr1;
begin
ha:=0.0;
for j:=1 to n do ha:=ha+H[ii,j]*A[j];
abc2:=ha;
// for ii:=1 to n do hdg[ii]:=abc2(ii,n,H,dg); 在DFPBFG程序中调用
end;
//===========================================================================
procedure sjxs_1;
var
ii : integer;
jg:array[1..400] of string;
begin
with form1.sumt do
begin
jg[1]:=floattostr(n);
jg[2]:=floattostr(kg);
jg[3]:=floattostr(kh);
jg[4]:=floattostr(R);
jg[5]:=floattostr(Cr);
jg[6]:=floattostr(T0);
jg[7]:=floattostr(EPS);
with form2.jgxs.lines do
begin
clear;
add(' 常用优化方法——惩罚函数法 ');
add(' ^^^^^^^^^^^^^^^^^^^^^^^^^^ ');
add(' ');
add('一、初始数据');
add('===============================================================================');
Add(' 设计变量个数 N = '+jg[1]);
Add(' -----------------------------------------------------------------------------');
Add(' 不等式约束个数 KG = '+jg[2]+' 等式约束个数 KH = '+jg[3]);
Add(' -----------------------------------------------------------------------------');
Add(' 惩罚因子 R = '+jg[4]+' 惩罚因子降低系数 C = '+jg[5]);
Add(' -----------------------------------------------------------------------------');
Add(' 初始步长 T0 = '+jg[6]+' 收敛精度 EPS = '+jg[7]);
Add(' -----------------------------------------------------------------------------');
Add(' 无约束优化方法: '+KCHECK_Z);
Add(' -----------------------------------------------------------------------------');
Add(' 设计变量初始点 X0:');
for ii := 1 to n do Add(#9+#9+#9+'X['+inttostr(ii)+']='+floattostr(x[ii]));
Add(' -----------------------------------------------------------------------------');
Add(' 设计变量下界 BL:');
for ii := 1 to n do Add(#9+#9+#9+'BL['+inttostr(ii)+']='+floattostr(BL[ii]));
Add(' -----------------------------------------------------------------------------');
Add(' 设计变量上界 BU:');
for ii := 1 to n do Add(#9+#9+#9+'BU['+inttostr(ii)+']='+floattostr(BU[ii]));
Add(' -----------------------------------------------------------------------------');
Add(' 初始点目标函数值 F(X0)= '+floattostr(FX));
if kg>0 then
begin
Add(' -----------------------------------------------------------------------------');
Add(' 初始点处的不等约束函数值 G(X0):');
for ii := 1 to kg do Add(#9+#9+#9+#9+'GX['+inttostr(ii)+']= '+formatfloat('0.000000E+00',GX[ii]));
end;
if kh>0 then
begin
Add(' -----------------------------------------------------------------------------');
Add(' 初始点处的等式约束函数值 H(X0):');
for ii := 1 to kh do Add(#9+#9+#9+#9+'HX['+inttostr(ii)+']= '+formatfloat('0.000000E+00',HX[ii]));
end;
end;
end;
end;
procedure sjxs_2;
var
ii : integer;
jg:array[1..400] of string;
begin
with form1.sumt do
begin
jg[1]:=inttostr(IRC);
jg[2]:=inttostr(ITE);
jg[3]:=inttostr(ILI);
jg[4]:=inttostr(NPE);
jg[5]:=inttostr(NFX);
with form2.jgxs.lines do
begin
add(' ');
add('三、优化结果__数据');
add('===============================================================================');
Add(' 罚函数构造次数 IRC = '+jg[1]);
Add(' -----------------------------------------------------------------------------');
Add(' 无约束优化方法调用次数 ITE = '+jg[2]+' 一维搜索方法调用次数 ILI = '+jg[3]);
Add(' -----------------------------------------------------------------------------');
Add(' 惩罚函数值计算次数 NPE = '+jg[4]+' 目标函数值计算次数 IFX = '+jg[5]);
Add(' -----------------------------------------------------------------------------');
Add(' 设计变量最优点 X*:');
for ii := 1 to n do Add(#9+#9+#9+'X['+inttostr(ii)+']= '+formatfloat('0.000000E+00',X[ii]));
Add(' -----------------------------------------------------------------------------');
Add(' 最优值 F(X*)= '+floattostr(FX));
if kg>0 then
begin
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;
if kh>0 then
begin
Add(' -----------------------------------------------------------------------------');
Add(' 最优点处的等式约束函数值 H(X*):');
for ii := 1 to kh do Add(#9+#9+#9+#9+'HX['+inttostr(ii)+']= '+formatfloat('0.000000E+00',HX[ii]));
end;
add('-------------------------------------------------------------------------------');
Add('--- STOP --- ');
end;
end;
end;
end.
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -