solutions for every iteration

2 views (last 30 days)
sura naji
sura naji on 23 Jan 2022
Answered: Image Analyst on 23 Jan 2022
format long
clear;
clc;
fun=@testcost
lb=[0.1,0.0003];
ub=[1.0,0.01];
%lb=[ 1,1,1,1,1,1]; %lower bound of variables
%ub=[10,10,10,10,10,10]; % upper bound of variables
nvars=2 % number rof variables
options = optimoptions('particleswarm')
options.SwarmSize=2000
options.Display='iter'
options = optimoptions(options,'PlotFcn',@pswplotbestf);
[x,fval,exitflag,output,MaxIterations] = particleswarm(fun,nvars,lb,ub,options);
%As=pi*(x(1)/2)^2-pi*((x(1)-2*x(2))/2)^2;
%Ac=pi*((x(1)-2*x(2))/2)^2;
%h=3.0;
%csteel=0.75*7850;
%cconc=180;
%cost=As*h*csteel+Ac*h*cconc
%****************************************************************
Ninp=6000;
Minp=100;
%****************************************************************
clc
[penalty,Mcalc,cost,x,Nmax,lambda, slenderness] = testcostbest(x,Ninp,Minp)
%diameter=x(1)
%thickness=x(2)
h=3;
fc=35;
fy=360;
Es=200;
Ec=((fc+8)/10)^0.3*22;
%-----------------calculations for global buckling ---------------
Asb=(pi*(x(1)/2)^2-pi*((x(1)-2*x(2))/2)^2);
Acb=(pi*((x(1)-2*x(2))/2)^2);
Npl=(fy*Asb+fc*Acb)*1000;
moicb=(1/4*pi*((x(1)-2*x(2))/2)^4);
moisb=((1/4*pi*(x(1)/2)^4)-(1/4*pi*((x(1)-2*x(2))/2)^4));
k=0.7;
EIeff=Es*moisb+0.6*Ec*moicb;
Ncr=pi^2*EIeff/(k*h)^2*1000000;
lambda=sqrt(Npl/Ncr);
%--------------slenderness-----------------------------------------
slenderness=x(1)/x(2);
% objective function
function f = testcost(x,Ninp,Minp)
% input axşal load (kN) and bending moment(kN.m)
%***********************************************************
Ninp=6000;
Minp=100;
%**********************************************************
% costs of steel and concrete
csteel=0.75*7850;
cconc=180;
% cross-sectional areas (mm2), moment of inertias (mm4), section modulus (mm3)
As=pi*(x(1)/2)^2-pi*((x(1)-2*x(2))/2)^2;
Ac=pi*((x(1)-2*x(2))/2)^2;
%moic=1/4*pi*((x(1)-2*x(2))/2)^4;
%mois=(1/4*pi*(x(1)/2)^4)-(1/4*pi*((x(1)-2*x(2))/2)^4);
sm=pi/32*(x(1)-2*x(2))^3;
% material properties: yield streses (MPa), leastic modulus (GPA) and height(m)
fc=35;
fy=360;
Es=200;
Ec=((fc+8)/10)^0.3*22;
h=3.0;
% other variables and constants
a1=817.2;
a2=5528;
a3=128.3;
a4=1282;
b1=595.2;
b2=5924;
b3=275.5;
b4=2103;
c1=22.5;
c2=-6358;
c3=6.827;
c4=-2257;
d1=385.7;
d2=1961;
d3=108.4;
d4=213.9;
teta=As*fy/(Ac*fc);
Nmax=Ac*fc*(a1+a2*x(2)/x(1)+a3*fc/fy+a4*teta);
Mmax=sm*fc*(b1+b2*x(2)/x(1)+b3*fc/fy+b4*teta);
M0=sm*fc*(c1+c2*x(2)/x(1)+c3*fc/fy-c4*teta);
NMmax=Ac*fc*(d1+d2*x(2)/x(1)+d3*fc/fy+d4*teta);
% coefficients of the polynomial
k1=M0;
k2=((2*Nmax^4-4*Nmax^2*NMmax^2)*(Mmax-M0)-2*M0*NMmax^2)/(Nmax*NMmax*(2*NMmax^3-3*Nmax*NMmax^2+Nmax^3));
k3=((4*Nmax*NMmax^3-Nmax^4)*(Mmax-M0)+3*M0*NMmax^4)/(Nmax*NMmax^2*(2*NMmax^3-3*Nmax*NMmax^2+Nmax^3));
k4=((Nmax^2-2*Nmax*NMmax)*(Mmax-M0)-M0*NMmax^2)/(Nmax*NMmax^2*(2*NMmax^3-3*Nmax*NMmax^2+Nmax^3));
Mcalc=k1+k2*Ninp+k3*Ninp^2+k4*Ninp^4;
%
% ---------penalties-------------
penalty=0;
%--------penalty #1 ------------(d-2*t>0)--------
cc1=x(1)-2*x(2);
if cc1<0
penalty=abs(cc1);
end
%-------penalty #2 -----------
cc2=x(1)/x(2);
if cc2>90*235/fy
penalty=penalty+100*abs(cc2);
end
slenderness=x(1)/x(2);
%--------penalty #3 ----------
if Minp>Mcalc
cc3=abs(Minp-Mcalc);
penalty=penalty+cc3*10;
end
%--------penalty #4 ----------
Asb=(pi*(x(1)/2)^2-pi*((x(1)-2*x(2))/2)^2);
Acb=(pi*((x(1)-2*x(2))/2)^2);
Npl=(fy*Asb+fc*Acb)*1000;
moicb=(1/4*pi*((x(1)-2*x(2))/2)^4);
moisb=((1/4*pi*(x(1)/2)^4)-(1/4*pi*((x(1)-2*x(2))/2)^4));
k=0.7;
EIeff=Es*moisb+0.6*Ec*moicb;
Ncr=pi^2*EIeff/(k*h)^2*1000000;
lambda=sqrt(Npl/Ncr);
cc4=lambda-2;
if cc4>0
penalty=penalty+abs(cc4)*1000;
end
%--------------------------------
% objective function
cost=As*h*csteel+Ac*h*cconc;
%cost1=abs(Minp-Mcalc)
f=cost*(1+penalty*100000);
if As<0
f=1000000000;
end
if Nmax<Ninp
f=1000000000;
penalty=f;
end
penalty;
end
% objective function
function [penalty,Mcalc,cost,x,Nmax, lambda, slenderness] = testcostbest(x, Ninp, Minp)
% costs of steel and concrete
csteel=0.75*7850;
cconc=180;
% cross-sectional areas (mm2), moment of inertias (mm4), section modulus (mm3)
As=pi*(x(1)/2)^2-pi*((x(1)-2*x(2))/2)^2;
Ac=pi*((x(1)-2*x(2))/2)^2;
%moic=1/4*pi*((x(1)-2*x(2))/2)^4;
%mois=(1/4*pi*(x(1)/2)^4)-(1/4*pi*((x(1)-2*x(2))/2)^4);
sm=pi/32*(x(1)-2*x(2))^3;
% material properties: yield streses (MPa), leastic modulus (GPA) and height(m)
fc=35;
fy=360;
Es=200;
Ec=((fc+8)/10)^0.3*22;
h=3.0;
% other variables and constants
a1=817.2;
a2=5528;
a3=128.3;
a4=1282;
b1=595.2;
b2=5924;
b3=275.5;
b4=2103;
c1=22.5;
c2=-6358;
c3=6.827;
c4=-2257;
d1=385.7;
d2=1961;
d3=108.4;
d4=213.9;
teta=As*fy/(Ac*fc);
Nmax=Ac*fc*(a1+a2*x(2)/x(1)+a3*fc/fy+a4*teta);
Mmax=sm*fc*(b1+b2*x(2)/x(1)+b3*fc/fy+b4*teta);
M0=sm*fc*(c1+c2*x(2)/x(1)+c3*fc/fy-c4*teta);
NMmax=Ac*fc*(d1+d2*x(2)/x(1)+d3*fc/fy+d4*teta);
% coefficients of the polynomial
k1=M0;
k2=((2*Nmax^4-4*Nmax^2*NMmax^2)*(Mmax-M0)-2*M0*NMmax^2)/(Nmax*NMmax*(2*NMmax^3-3*Nmax*NMmax^2+Nmax^3));
k3=((4*Nmax*NMmax^3-Nmax^4)*(Mmax-M0)+3*M0*NMmax^4)/(Nmax*NMmax^2*(2*NMmax^3-3*Nmax*NMmax^2+Nmax^3));
k4=((Nmax^2-2*Nmax*NMmax)*(Mmax-M0)-M0*NMmax^2)/(Nmax*NMmax^2*(2*NMmax^3-3*Nmax*NMmax^2+Nmax^3));
Mcalc=k1+k2*Ninp+k3*Ninp^2+k4*Ninp^4;
%
% ---------penalties-------------
penalty=0;
%--------penalty #1 ------------(d-2*t>0)--------
cc1=x(1)-2*x(2);
if cc1<0
penalty=abs(cc1);
end
%-------penalty #2 -----------
cc2=x(1)/x(2);
if cc2>90*235/fy
penalty=penalty+100*abs(cc2);
end
slenderness=x(1)/x(2);
%--------penalty #3 ----------
if Minp>Mcalc
cc3=abs(Minp-Mcalc);
penalty=penalty+cc3*10;
end
%--------penalty #4 ----------
Asb=(pi*(x(1)/2)^2-pi*((x(1)-2*x(2))/2)^2);
Acb=(pi*((x(1)-2*x(2))/2)^2);
Npl=(fy*Asb+fc*Acb)*1000;
moicb=(1/4*pi*((x(1)-2*x(2))/2)^4);
moisb=((1/4*pi*(x(1)/2)^4)-(1/4*pi*((x(1)-2*x(2))/2)^4));
k=0.7;
EIeff=Es*moisb+0.6*Ec*moicb;
Ncr=pi^2*EIeff/(k*h)^2*1000000;
lambda=sqrt(Npl/Ncr);
cc4=lambda-2;
if cc4>0
penalty=penalty+abs(cc4)*1000;
end
%--------------------------------
% objective function
cost=As*h*csteel+Ac*h*cconc;
%cost1=abs(Minp-Mcalc)
f=cost*(1+penalty*100000);
if As<0
f=1000000000;
end
if Nmax<Ninp
f=1000000000;
penalty=f;
end
penalty;
end
Can anyone help me please, when l do run l get just the optimum solutions, l need the cost and valued x1 and x2 for evey iteration ?

Sign in to answer this question.

Answers (1)

Image Analyst
Image Analyst on 23 Jan 2022
You need to index your variables to make them vectors (lists of values) rather than just overwriting them every iterations, like instead of
myVariable = whatever;
do
myVariable(index) = whatever;
where index is your "for" or "while" loop iterator variable that gets incremented on each iteration.

Categories

Find more on Whos in Help Center and File Exchange

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!