You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
311 lines
10 KiB
311 lines
10 KiB
% concurrent:
|
|
% loop = 83;
|
|
% compliance=923.53;
|
|
%
|
|
% c-a-b:
|
|
% loop=101;
|
|
% compliance = 923.40;
|
|
|
|
function ani_opt_concurrent_smallrange()
|
|
nx=60;
|
|
ny=1;
|
|
nz=20;
|
|
volfrac=0.5;
|
|
rmin=1.5;
|
|
[Hs,H]=filter_coef_3D(nx,ny,nz,rmin);
|
|
|
|
a(1:nz,1:ny,1:nx) = 1.0;
|
|
b(1:nz,1:ny,1:nx) = 1.0;
|
|
c(1:nz,1:ny,1:nx) = 0.0;
|
|
aPhys = a;
|
|
bPhys = b;
|
|
cPhys = c;
|
|
% -1+a+b-c>=0
|
|
% 1-a+b-c>=0
|
|
% 1+a-b-c>=0
|
|
|
|
%OPTIMIZE a
|
|
% a>=max(1-b+c,-1+b+c)
|
|
% a<=1+b-c
|
|
%OPTIMIZE b
|
|
% b>=max(1-a+c,-1+a+c)
|
|
% b<=1+a-c
|
|
%OPTIMIZE c
|
|
% c<=min(-1+a+b,1-a+b,1+a-b)
|
|
|
|
m = 1; % The number of general constraints.
|
|
n = 3 * nx * ny * nz; % The number of design variables x_j.
|
|
nele = nx * ny * nz;
|
|
xmin = [ones(nele,1).*0.6;...
|
|
ones(nele,1).*0.6;...
|
|
ones(nele,1).*(-4.9)];
|
|
% Column vector with the lower bounds for the variables x_j.
|
|
xmax = [ones(nele,1).*1.9;...
|
|
ones(nele,1).*1.9;...
|
|
ones(nele,1).*0.9];
|
|
xold1 = [ a(:); b(:);c(:)]; % xval, one iteration ago (provided that iter>1).
|
|
xold2 = [ a(:); b(:);c(:)]; % xval, two iterations ago (provided that iter>2).
|
|
low = ones(n,1); % Column vector with the lower asymptotes from the previous iteration (provided that iter>1).
|
|
upp = ones(n,1); % Column vector with the upper asymptotes from the previous iteration (provided that iter>1).
|
|
a0 = 1; % The constants a_0 in the term a_0*z.
|
|
a1 = zeros(m,1); % Column vector with the constants a_i in the terms a_i*z.
|
|
x_MMA = 3e3.*ones(m,1); % Column vector with the constants c_i in the terms c_i*y_i.
|
|
d = zeros(m,1); % Column vector with the constants d_i in the terms 0.5*d_i*(y_i)^2.
|
|
|
|
|
|
maxloop = 300;
|
|
a_log = cell(maxloop,2);
|
|
b_log = cell(maxloop,2);
|
|
c_log = cell(maxloop,2);
|
|
tim=zeros(maxloop,1);
|
|
sumCompliance=zeros(maxloop,1);
|
|
volfraction = zeros(maxloop,1);
|
|
loop = 0;
|
|
change = 1.;
|
|
while change > 0.001 && loop < maxloop
|
|
loop = loop + 1;
|
|
tic
|
|
|
|
[a,b,c,aPhys,bPhys,cPhys,xold1,xold2,change,compliance,vol,low,upp,comp]=...
|
|
optimizeABC(a,b,c,aPhys,bPhys,cPhys,xold1,xold2,Hs,H,volfrac,low,upp,nx,ny,nz,loop,a0,a1,x_MMA,d,...
|
|
n,m,xmin,xmax);
|
|
sumCompliance(loop,1)=compliance;
|
|
volfraction(loop,1) = mean(vol);
|
|
disp([' It.: ' sprintf('%4i',loop)...
|
|
' compliance.: ' sprintf('%10.10f',compliance)...
|
|
' volfrac:' sprintf('%10.10f',mean(vol))...
|
|
' change.: ' sprintf('%6.3f',change )]);
|
|
a_log{loop,1}=a;a_log{loop,2}=aPhys;
|
|
b_log{loop,1}=b;b_log{loop,2}=bPhys;
|
|
c_log{loop,1}=c;c_log{loop,2}=cPhys;
|
|
|
|
% plot_comp(comp);
|
|
|
|
% plot_distributions(aPhys,bPhys,cPhys);
|
|
plot_convergence_all(loop,maxloop,sumCompliance,volfraction);
|
|
|
|
tim(loop)=toc;
|
|
save('ani_rmin1.5_concurrent_oldH.mat','aPhys','bPhys','cPhys',...
|
|
'sumCompliance','volfraction','a_log','b_log','c_log','tim');
|
|
% save ani_DATA_cab.mat
|
|
end
|
|
disp(['Time = ',num2str(sum(tim))])
|
|
% ch=interp_ch(aPhys,bPhys,cPhys);
|
|
% [ch_dc]=interp_ch_dc(aPhys,bPhys,cPhys);
|
|
% ke_all=lk_all(ch);ke_dc=lk_all(ch_dc);
|
|
% U=FEM_Beam(nx,ny,nz,ke_all);
|
|
% [comp,compliance,dcdc] = COMP(U,ke_all,ke_dc,nx,ny,nz);
|
|
|
|
disp(['C = ',num2str(compliance)])
|
|
figure()
|
|
set(gcf,'position',[100 650 300 300])
|
|
cc = squeeze(comp(:,1,:));
|
|
colormap(1-gray); imagesc(cc);
|
|
axis equal; axis tight; axis off;
|
|
% title('Compliance Distribution')
|
|
set(gca,'YDir','normal')
|
|
end
|
|
|
|
function plot_convergence_all(loop,maxloop,sumCompliance,volfraction)
|
|
figure(3)
|
|
% figure()
|
|
% loop=75;maxloop=loop;
|
|
set(gcf,'position',[100 100 1200 400])
|
|
subplot(2,1,1)
|
|
plotC=sumCompliance.';
|
|
plot([1:loop].',plotC(1:loop))
|
|
xlim([0 maxloop])
|
|
subplot(2,1,2)
|
|
plotV=volfraction.';
|
|
plot([1:loop].',plotV(1:loop))
|
|
xlim([0 maxloop])
|
|
pause(1e-6);
|
|
end
|
|
|
|
function [a,b,c,aPhys,bPhys,cPhys,xold1,xold2,change,compliance,vol,low,upp,comp]=...
|
|
optimizeABC(a,b,c,aPhys,bPhys,cPhys,xold1,xold2,Hs,H,volfrac,low,upp,nx,ny,nz,loop,a0,a1,x_MMA,d,...
|
|
n,m,xmin,xmax)
|
|
nele = nx * ny * nz;
|
|
ch=interp_ch(aPhys,bPhys,cPhys);
|
|
ke_all=lk_all(ch);
|
|
U=FEM_Beam(nx,ny,nz,ke_all);
|
|
|
|
[ch_da]=interp_ch_da(aPhys,bPhys,cPhys);
|
|
[ch_db]=interp_ch_db(aPhys,bPhys,cPhys);
|
|
[ch_dc]=interp_ch_dc(aPhys,bPhys,cPhys);
|
|
ke_da=lk_all(ch_da);
|
|
ke_db=lk_all(ch_db);
|
|
ke_dc=lk_all(ch_dc);
|
|
[~,~,dcdb] = COMP(U,ke_all,ke_db,nx,ny,nz);
|
|
[~,~,dcda] = COMP(U,ke_all,ke_da,nx,ny,nz);
|
|
[comp,compliance,dcdc] = COMP(U,ke_all,ke_dc,nx,ny,nz);
|
|
|
|
[dvda]=interp_dvda(aPhys(:),bPhys(:),cPhys(:));
|
|
[dvdb]=interp_dvdb(aPhys(:),bPhys(:),cPhys(:));
|
|
[dvdc]=interp_dvdc(aPhys(:),bPhys(:),cPhys(:));
|
|
|
|
dv=[dvda;dvdb;dvdc];
|
|
xval=[aPhys(:);bPhys(:);cPhys(:)];
|
|
f0val = compliance;
|
|
df0dx = [dcda(:);dcdb(:);dcdc(:)];
|
|
vol=interp_v(aPhys(:),bPhys(:),cPhys(:));
|
|
fval = sum(vol)/(volfrac*nele)-1;
|
|
dfdx = dv(:)' / (volfrac*nele);
|
|
|
|
|
|
[xmma, ~, ~, ~, ~, ~, ~, ~, ~, low,upp] = ...
|
|
mmasub(m, n, loop, xval, xmin, xmax, xold1, xold2, ...
|
|
f0val,df0dx,fval,dfdx,low,upp,a0,a1,x_MMA,d);
|
|
|
|
a = reshape(xmma(1:nele),nz,ny,nx);
|
|
b = reshape(xmma((nele+1):2*nele),nz,ny,nx);
|
|
c = reshape(xmma((2*nele+1):n),nz,ny,nx);
|
|
aPhys(:) = (double(H)*double(a(:)))./double(Hs);
|
|
bPhys(:) = (double(H)*double(b(:)))./double(Hs);
|
|
cPhys(:) = (double(H)*double(c(:)))./double(Hs);
|
|
xold2 = xold1(:);
|
|
xold1 = xval(:);
|
|
|
|
change = max(max(abs([a(:);b(:);c(:)]-xval)));
|
|
|
|
end
|
|
|
|
function [c,cPhys,cold1,cold2,change,compliance,vol,clow,cupp,comp]=...
|
|
optimizeC(a,b,c,aPhys,bPhys,cPhys,cold1,cold2,Hs,H,volfrac,clow,cupp,nx,ny,nz,loop)
|
|
m = 1; % The number of general constraints.
|
|
n = nx * ny * nz; % The number of design variables x_j.
|
|
nele = nx * ny * nz;
|
|
a0 = 1; % The constants a_0 in the term a_0*z.
|
|
a1 = zeros(m,1); % Column vector with the constants a_i in the terms a_i*z.
|
|
x_MMA = 9.5e2.*ones(m,1); % Column vector with the constants c_i in the terms c_i*y_i.
|
|
d = zeros(m,1); % Column vector with the constants d_i in the terms 0.5*d_i*(y_i)^2.
|
|
|
|
ch=interp_ch(aPhys,bPhys,cPhys);
|
|
ke_all=lk_all(ch);
|
|
U=FEM_Beam(nx,ny,nz,ke_all);
|
|
|
|
%OPTIMIZE c
|
|
% c<=min(-1+a+b,1-a+b,1+a-b)
|
|
% c_min = max(1-a+c,-1+a+c);
|
|
% c_min = -1-a-b;
|
|
% c_max = min(min(-1+a+b,1-a+b),1+a-b);
|
|
c_min = -1-aPhys-bPhys;
|
|
c_max = min(min(-1+aPhys+bPhys,1-aPhys+bPhys),1+aPhys-bPhys);
|
|
xmin = max(c_min(:),-4.98);% Column vector with the lower bounds for the variables x_j.
|
|
xmax = min(c_max(:),0.98);
|
|
|
|
[ch_dc]=interp_ch_dc(aPhys,bPhys,cPhys);
|
|
ke_dc=lk_all(ch_dc);
|
|
[comp,compliance,dcdc] = COMP(U,ke_all,ke_dc,nx,ny,nz);
|
|
[dvdc]=interp_dvdc(aPhys(:),bPhys(:),cPhys(:));
|
|
dv=dvdc;
|
|
xval=cPhys(:);
|
|
f0val = compliance;
|
|
df0dx = dcdc(:);
|
|
vol=interp_v(aPhys(:),bPhys(:),cPhys(:));
|
|
fval = sum(vol)/(volfrac*nele)-1;
|
|
dfdx = dv(:)' / (volfrac*nele);
|
|
[xmma, ~, ~, ~, ~, ~, ~, ~, ~, clow,cupp] = ...
|
|
mmasub(m, n, loop, xval, xmin, xmax, cold1, cold2, ...
|
|
f0val,df0dx,fval,dfdx,clow,cupp,a0,a1,x_MMA,d);
|
|
c = reshape(xmma,nz,ny,nx);
|
|
cPhys(:) = (double(H)*double(c(:)))./double(Hs);
|
|
cold2 = cold1(:);
|
|
cold1 = xval(:);
|
|
change = max(max(abs(c(:)-xval)));
|
|
vol=interp_v(aPhys(:),bPhys(:),cPhys(:));
|
|
end
|
|
|
|
|
|
|
|
|
|
function [b,bPhys,bold1,bold2,change,compliance,vol,blow,bupp,comp]=...
|
|
optimizeB(a,b,c,aPhys,bPhys,cPhys,bold1,bold2,Hs,H,volfrac,blow,bupp,nx,ny,nz,loop)
|
|
m = 1; % The number of general constraints.
|
|
n = nx * ny * nz; % The number of design variables x_j.
|
|
nele = nx * ny * nz;
|
|
a0 = 1; % The constants a_0 in the term a_0*z.
|
|
a1 = zeros(m,1); % Column vector with the constants a_i in the terms a_i*z.
|
|
x_MMA = 1e2.*ones(m,1); % Column vector with the constants c_i in the terms c_i*y_i.
|
|
d = zeros(m,1); % Column vector with the constants d_i in the terms 0.5*d_i*(y_i)^2.
|
|
|
|
ch=interp_ch(aPhys,bPhys,cPhys);
|
|
ke_all=lk_all(ch);
|
|
U=FEM_Beam(nx,ny,nz,ke_all);
|
|
|
|
%OPTIMIZE b
|
|
% b>=max(1-a+c,-1+a+c)
|
|
% a<=1+a-c
|
|
% b_min = max(max(1-a+c,-1+a+c),-1-a-c);
|
|
% b_max = 1+a-c;
|
|
b_min = max(max(1-aPhys+cPhys,-1+aPhys+cPhys),-1-aPhys-cPhys);
|
|
b_max = 1+aPhys-cPhys;
|
|
xmin = max(b_min(:),0.52);% Column vector with the lower bounds for the variables x_j.
|
|
xmax = min(b_max(:),1.98);
|
|
|
|
[ch_db]=interp_ch_db(aPhys,bPhys,cPhys);
|
|
ke_db=lk_all(ch_db);
|
|
[comp,compliance,dcdb] = COMP(U,ke_all,ke_db,nx,ny,nz);
|
|
[dvdb]=interp_dvdb(aPhys(:),bPhys(:),cPhys(:));
|
|
dv=dvdb;
|
|
xval=bPhys(:);
|
|
f0val = compliance;
|
|
df0dx = dcdb(:);
|
|
vol=interp_v(aPhys(:),bPhys(:),cPhys(:));
|
|
fval = sum(vol)/(volfrac*nele)-1;
|
|
dfdx = dv(:)' / (volfrac*nele);
|
|
[xmma, ~, ~, ~, ~, ~, ~, ~, ~, blow,bupp] = ...
|
|
mmasub(m, n, loop, xval, xmin, xmax, bold1, bold2, ...
|
|
f0val,df0dx,fval,dfdx,blow,bupp,a0,a1,x_MMA,d);
|
|
b = reshape(xmma,nz,ny,nx);
|
|
bPhys(:) = (double(H)*double(b(:)))./double(Hs);
|
|
bold2 = bold1(:);
|
|
bold1 = xval(:);
|
|
change = max(max(abs(b(:)-xval)));
|
|
vol=interp_v(aPhys(:),bPhys(:),cPhys(:));
|
|
end
|
|
|
|
function [a,aPhys,aold1,aold2,change,compliance,vol,alow,aupp,comp]=...
|
|
optimizeA(a,b,c,aPhys,bPhys,cPhys,aold1,aold2,Hs,H,volfrac,alow,aupp,nx,ny,nz,loop)
|
|
m = 1; % The number of general constraints.
|
|
n = nx * ny * nz; % The number of design variables x_j.
|
|
nele = nx * ny * nz;
|
|
|
|
a0 = 1; % The constants a_0 in the term a_0*z.
|
|
a1 = zeros(m,1); % Column vector with the constants a_i in the terms a_i*z.
|
|
x_MMA = 1e2.*ones(m,1); % Column vector with the constants c_i in the terms c_i*y_i.
|
|
d = zeros(m,1); % Column vector with the constants d_i in the terms 0.5*d_i*(y_i)^2.
|
|
|
|
ch=interp_ch(aPhys,bPhys,cPhys);
|
|
ke_all=lk_all(ch);
|
|
U=FEM_Beam(nx,ny,nz,ke_all);
|
|
|
|
% a_min = max(max(1-b+c,-1+b+c),-1-b-c);
|
|
% a_max = 1+b-c;
|
|
a_min = max(max(1-bPhys+cPhys,-1+bPhys+cPhys),-1-bPhys-cPhys);
|
|
a_max = 1+bPhys-cPhys;
|
|
xmin = max(a_min(:),0.52);% Column vector with the lower bounds for the variables x_j.
|
|
xmax = min(a_max(:),1.98);
|
|
|
|
[ch_da]=interp_ch_da(aPhys,bPhys,cPhys);
|
|
ke_da=lk_all(ch_da);
|
|
[comp,compliance,dcda] = COMP(U,ke_all,ke_da,nx,ny,nz);
|
|
|
|
[dvda]=interp_dvda(aPhys(:),bPhys(:),cPhys(:));
|
|
dv=dvda;
|
|
xval=aPhys(:);
|
|
f0val = compliance;
|
|
df0dx = dcda(:);
|
|
vol=interp_v(aPhys(:),bPhys(:),cPhys(:));
|
|
fval = sum(vol)/(volfrac*nele)-1;
|
|
dfdx = dv(:)' / (volfrac*nele);
|
|
[xmma, ~, ~, ~, ~, ~, ~, ~, ~, alow,aupp] = ...
|
|
mmasub(m, n, loop, xval, xmin, xmax, aold1, aold2, ...
|
|
f0val,df0dx,fval,dfdx,alow,aupp,a0,a1,x_MMA,d);
|
|
a = reshape(xmma,nz,ny,nx);
|
|
aPhys(:) = (double(H)*double(a(:)))./double(Hs);
|
|
aold2 = aold1(:);
|
|
aold1 = xval(:);
|
|
change = max(max(abs(a(:)-xval)));
|
|
vol=interp_v(aPhys(:),bPhys(:),cPhys(:));
|
|
end
|