📄 zhj.m
字号:
% This procedure can optimize an H-type section steel pole ,
% which is worked out by Zhang Haojie ( 0302153 )
% 2004/12/11
% Units: N,mm,Mpa
% main function : SUMT------------------------------------------
function opt = optimum(X0,r,ep1,ep2,C)
global P Segma l E a1 a2 mu
global N
global ep1 ep2
% input data initialization
P=1.5e6;
Segma=350;
l=2e4;
%l=3e3;
E=2.1e5;
a1=0.5;
a2=0.5;
mu=1;
%
X0=[1,1,40,500]; % l=2e4 r=[1e6,1e6,1e6,1e6,1e6,1e4,1e4,1e5,1e6] h=1e-3; r=1e7*[1,1,1,1,1,1,1,1,1] h=1e-3
% l=3e3 r=[1e6,1e7,1e7,1e7,1e7,1e4,1e4,1e5,1e6] h=1e-3 ; r=1e7*[1,1,1,1,1,1,1,1,1] h=1e-6
%X0=[2,3,15,300]; % l=2e4 r=[1e6,1e6,1e6,1e6,1e6,1e4,1e4,1e5,1e6] h=1e-3 ; r=1e7*[1,1,1,1,1,1,1,1,1] h=1e-3
% l=3e3 r=[1e6,1e7,1e7,1e6,1e6,1e4,1e4,1e5,1e6] h=1e-3 ; r=1e7*[1,1,1,1,1,1,1,1,1] h=1e-3
%X0=[2,2,10,400]; % l=2e4 r=[1e6,1e6,1e6,1e6,1e6,1e4,1e4,1e5,1e6] h=1e-3
% l=3e3 r=[1e6,1e7,1e7,1e6,1e6,1e4,1e4,1e5,1e6] h=1e-5
%X0=[2,3,20,500]; % l=2e4 r=[1e6,1e6,1e6,1e6,1e6,1e4,1e4,1e5,1e6] h=1e-3
% l=3e3 r=[1e6,1e8,1e8,1e6,1e6,1e4,1e4,1e5,1e6] h=1e-5
%X0=[2,3,30,500]; % l=2e4 r=[1e6,1e7,1e7,1e7,1e7,1e4,1e4,1e5,1e6] h=1e-3
% l=3e3 r=[1e6,1e8,1e8,1e8,1e8,1e4,1e4,1e5,1e6] h=1e-4
X1=X0;
X2=X0;
N=4; %dimension of the problem
r=[1e6,1e6,1e6,1e6,1e6,1e4,1e4,1e5,1e6];
%r=1e7*[1,1,1,1,1,1,1,1,1];
ep1=1e-3;
ep2=1e-2;
C=0.9;
k=0;
flag=1;
% start to estimate elapsed time
tic
% SUMT procedure
while flag==1
k=k+1;
X2=Powell(X1,r);
flag=judge(X1,X2);
if flag==1
r=C*r;
X1=X2;
end
end
% output data
k
X2
opt=A(X2)
seg=P/A(X2)
gu1=seg-Segma
gu2=seg-(pi^2*E/(mu*l)^2)*(X2(1)^3*X2(2)*X2(4)^2*a1/(6*(1+2*X2(1)*X2(2))))
gu3=seg-(pi^2*E/(mu*l)^2)*((1+6*X2(1)*X2(2))*X2(4)^2*a1/(12*(1+2*X2(1)*X2(2))))
gu4=seg-3.62*E*(X2(3)/X2(4))^2*a2
gu5=seg-0.385*E*(2*X2(2)*X2(3)/(X2(1)*X2(4)))^2*a2
% time out
toc
% sub function
% object & constraints------------------------------------------
% object
function f=A(X)
f=(1+2*X(1)*X(2))*X(3)*X(4);
% constraint 1
function f=g1(X)
global P Segma l E a1 a2 mu
f=P/A(X)-Segma;
% constraint 2
function f=g2(X)
global P Segma l E a1 a2 mu
f=P/A(X)-(pi^2*E/(mu*l)^2)*(X(1)^3*X(2)*X(4)^2*a1/(6*(1+2*X(1)*X(2))));
% constraint 3
function f=g3(X)
global P Segma l E a1 a2 mu
f=P/A(X)-(pi^2*E/(mu*l)^2)*((1+6*X(1)*X(2))*X(4)^2*a1/(12*(1+2*X(1)*X(2))));
% constraint 4
function f=g4(X)
global P Segma l E a1 a2 mu
f=P/A(X)-3.62*E*(X(3)/X(4))^2*a2;
% constraint 5
function f=g5(X)
global P Segma l E a1 a2 mu
f=P/A(X)-0.385*E*(2*X(2)*X(3)/(X(1)*X(4)))^2*a2;
% constraint 6
function f=g6(X)
f=-X(1)+0.2;
% constraint 7
function f=g7(X)
f=-X(2)+0.2;
% constraint 8
function f=g8(X)
f=-X(3)+2;
% constraint 9
function f=g9(X)
f=-X(4)+200;
% punishment function-------------------------------------------
function f=puf(X,r)
g = [ 1/g1(X),1/g2(X),1/g3(X),1/g4(X),1/g5(X),1/g6(X),1/g7(X),1/g8(X),1/g9(X) ]';
f = A(X)-r*g;
% judgement function--------------------------------------------
% judge
function f=judge(X1,X2,ep1,ep2)
global N
global ep1 ep2
asum=0;
for k=1:1:N
asum=asum+(X1(k)-X2(k))^2;
end
asum=sqrt(asum);
f1=A(X1);
f2=A(X2);
if ~(asum>ep1)&~(abs((f1-f2)/f2)>ep2)
f=0;
else
f=1;
end
% judge1
function f=judge1(X1,X2,r)
global N
global ep1 ep2
asum=0;
for k=1:1:N
asum=asum+(X1(k)-X2(k))^2;
end
f1=puf(X1,r);
f2=puf(X2,r);
if ~(asum>ep1)&~(abs((f1-f2)/f2)>ep2)
f=0;
else
f=1;
end
% judge2
function f=judge2(f1,f2,f3,deltm)
if (f3<f1)&((f1-2*f2+f3)*(f1-f2-deltm)^2 < 0.5*deltm*(f1-f3)^2)
f=1;
else
f=0;
end
% Powell--------------------------------------------------------
function X2=Powell(X1,r)
global N
X2=X1;
Xa=[0,0,0,0;0,0,0,0;0,0,0,0;0,0,0,0;0,0,0,0;0,0,0,0];
Xa(1,:)=X1;
Xb=X1;
X=[0,0,0,0];
S1=[1,0,0,0;0,1,0,0;0,0,1,0;0,0,0,1;0,0,0,0];
S2=S1;
n=0;
flag1=1;
while flag1==1
n=n+1;
for k=1:1:N
Xa(k+1,:)=ods(Xa(k,:),S1(k,:),r);
end
S1(N+1,:)=Xa(N+1,:)-Xa(1,:);
Xa(N+2,:)=2*Xa(N+1,:)-Xa(1,:);
m=Max(Xa,r);
deltm=puf(Xa(m,:),r)-puf(Xa(m+1,:),r);
Sm=Xa(m+1,:)-Xa(m,:);
%
f1=puf(Xa(1,:),r);
f2=puf(Xa(N+1,:),r);
f3=puf(Xa(N+2,:),r);
flag2=judge2(f1,f2,f3,deltm);
%
if flag2==1
X=ods(Xa(N+1,:),S1(N+1,:),r);
Xb=X;
for ii=1:1:N
if ii< m
S2(ii,:)=S1(ii,:);
else
S2(ii,:)=S1(ii+1,:);
end
end
else
if f2<f3
Xb=Xa(N+1,:);
else
Xb=Xa(N+2,:);
end
S2=S1;
end
flag1=judge1(Xa(1,:),Xb,r);
Xa(1,:)=Xb;
S1=S2;
end
X2=Xa(1,:);
% deltm position -----------------------------------------------
function m=Max(Xa,r)
global N
m=1;
deltm=puf(Xa(1,:),r)-puf(Xa(2,:),r);
for k=2:1:N
delt=puf(Xa(k,:),r)-puf(Xa(k+1,:),r);
if deltm<delt
m=k;
end
end
% one dimension search function---------------------------------
function X2=ods(X1,S,r)
% determine the search space
alf0=0;
h=1e-3;
flag1=1;
flag2=1;
n=0;
alf1=alf0;
alf2=alf0+h;
Xa=X1+alf1*S;
Xb=X1+alf2*S;
f1=puf(Xa,r);
f2=puf(Xb,r);
if f2<f1
while flag1==1
n=n+1;
h=2*h;
alf2=alf2+h;
Xa=Xb;
Xb=X1+alf2*S;
f1=f2;
f2=puf(Xb,r);
if f2<f1
alf1=alf2-h;
else
a=alf1;
b=alf2;
flag1=0;
end
end
else
h=-h/4;
while flag2==1
n=n+1;
alf1=alf1+h;
Xb=Xa;
Xa=X1+alf1*S;
f2=f1;
f1=puf(Xa,r);
if f2<f1
a=alf1;
b=alf2;
flag2=0;
else
alf2=alf1-h;
h=2*h;
end
end
end
%determine the optimum search stepsize
ep=1e-4;
lda=0.618;
k=log(ep/(b-a))/log(lda)+1;
k=ceil(k);
alf1=b-lda*(b-a);
alf2=a+lda*(b-a);
Xa=X1+alf1*S;
Xb=X1+alf2*S;
f1=puf(Xa,r);
f2=puf(Xb,r);
flag=1;
while flag==1
if f1>f2
a=alf1;
alf1=alf2;
alf2=a+lda*(b-a);
Xa=Xb;
Xb=X1+alf2*S;
f1=f2;
f2=puf(Xb,r);
else
b=alf2;
alf2=alf1;
alf1=b-lda*(b-a);
Xb=Xa;
Xa=X1+alf1*S;
f2=f1;
f1=puf(Xa,r);
end
if k>1
k=k-1;
else
flag=0;
end
end
if f1<f2
alf=alf1;
else
alf=alf2;
end
X2=X1+alf*S;
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -