How to do optimization coding for minimize peak width ?

2024 3 20
1 个回答
回答已采纳
更新时间：2024 3 26
7 次查看（30 天）
Ran in:

0 个投票

Ran in:
i need a optimization code for minimize the "peak_width" (width at half maximum) by the optimized "del" (thickness variation parameter) that varying between -1 and 1 at specific frequency and i need unit transmission(maximum transmission) also (red solid line)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This is a program for 1-D Photonic crystal
% For Photonic band structure calculation and Transmission Spectrum
% This program uses Transfer Matrix Method
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc
clear
close all
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define material properties and simulation parameters
P = 0; %%% Hydrostatic pressure in MPa
n0 = 1; % Refractive index of air
ns = 1.46; % Refractive index of subtrate
n01 = 1.578; % Refractive index of PS
n02 = 1.484; % Refractive index of PMMA
% Calculate Epsilon for n01, n02
e01 = n01^2;
e02 = n02^2;
% Constants
p11_1 = 0.32;
p12_1 = 0.31;
p11_2 = 0.3;
p12_2 = 0.297;
v1 = 0.35;
E1 = 3.3e3; % Convert to scientific notation (3.3 x 10^3)
v2 = 0.37;
E2 = 3.0303e3; % Convert to scientific notation (3.0303 x 10^3)
% Formula for calculating nA
e1 = e01 - ((e01^2)/2) * (((p11_1/E1) * (v1+1) * P) + ((p12_1/E1) * (3*v1+1) * P));
nA = sqrt(e1);
% Formula for calculating nA
e2 = e02 - ((e02^2)/2) * (((p11_2/E2) * (v2+1) * P) + ((p12_2/E2) * (3*v2+1) * P));
nB = sqrt(e2);
del=0.131; %%% Thickness variation parameter
% Initialization of parameters
dair=1200e-9; % Thickness of air in meters
dA=780e-9; % Thickness of First Layer in meters
dB=830e-9; % Thickness of Second Layer in meters
ddA=(1+del)*dA; % Thickness of First Layer with thickness variation in meters
ddB=(1+del)*dB; % Thickness of Second Layer with thickness variation in meters
c=3e8; % Velocity of Light (m/s)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Frequency Range
% f=linspace(58,67,16001); % Frequncy in THz
f=58:0.001:67; % Frequncy in THz
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for loop=1:length(f)
%%% angular frequency
w=2*pi*f*1e12;
a=dA+dB; %%% Lattice Constant
DA=((w(loop))/c)*dA*nA; %%% wave number of first layer
DB=((w(loop))/c)*dB*nB; %%% wave number of second layer
Da=((w(loop))/c)*(dA*0.5)*nA; %%% wave number of first half layer
Db=((w(loop))/c)*(dB*0.5)*nB; %%% wave number of second half layer
DDA=((w(loop))/c)*ddA*nA; %%% wave number of first layer with thickness variation
DDB=((w(loop))/c)*ddB*nB; %%% wave number of second layer with thickness variation
DDa=((w(loop))/c)*(ddA*0.5)*nA; %%% wave number of first half layer with thickness variation
DDb=((w(loop))/c)*(ddB*0.5)*nB; %%% wave number of second half layer with thickness variation
Dair=((w(loop))/c)*(dair*0.5)*n0; %%% wave number of air layer
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Transfer Matrix elements of air
a11=cos(Dair); a12=-1i*sin(Dair)/n0; a21=-1i*n0*sin(Dair); a22=cos(Dair);
AA=[a11 a12; a21 a22];
%%% Transfer Matrix elements of first layer
m11=cos(DA); m12=-1i*sin(DA)/nA; m21=-1i*nA*sin(DA); m22=cos(DA);
mA=[m11 m12; m21 m22];
%%% Transfer MAtrix elements of Second layer
l11=cos(DB); l12=-1i*sin(DB)/nB; l21=-1i*nB*sin(DB); l22=cos(DB);
mB=[l11 l12; l21 l22];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Transfer Matrix elements of first half layer
ma11=cos(Da); ma12=-1i*sin(Da)/nA; ma21=-1i*nA*sin(Da); ma22=cos(Da);
ma=[ma11 ma12; ma21 ma22];
%%% Transfer MAtrix elements of Second half layer
lb11=cos(Db); lb12=-1i*sin(Db)/nB; lb21=-1i*nB*sin(Db); lb22=cos(Db);
mb=[lb11 lb12; lb21 lb22];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Transfer Matrix elements of first layer with thickness variation
dm11=cos(DDA); dm12=-1i*sin(DDA)/nA; dm21=-1i*nA*sin(DDA); dm22=cos(DDA);
mA1=[dm11 dm12; dm21 dm22];
%%% Transfer MAtrix elements of Second layer with thickness variation
dl11=cos(DDB); dl12=-1i*sin(DDB)/nB; dl21=-1i*nB*sin(DDB); dl22=cos(DDB);
mB2=[dl11 dl12; dl21 dl22];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Transfer Matrix elements of first half layer with thickness variation
dma11=cos(DDa); dma12=-1i*sin(DDa)/nA; dma21=-1i*nA*sin(DDa); dma22=cos(DDa);
ma1=[dma11 dma12; dma21 dma22];
%%% Transfer MAtrix elements of Second half layer with thickness variation
dlb11=cos(DDb); dlb12=-1i*sin(DDb)/nB; dlb21=-1i*nB*sin(DDb); dlb22=cos(DDb);
mb2=[dlb11 dlb12; dlb21 dlb22];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Unit cell
m=(mA*mB);
M0=(m^25); %%% A-B structure
M1=(mb*mA*mb)^25; %%% B/2-A-B/2 structure
M2=(ma*mB*ma)^25; %%% A/2-B-A/2 structure
M3=(mb2*mA1*mb2)^25; %%% B/2-A-B/2 structure with thickness variation
M4=(ma1*mB2*ma1)^25; %%% A/2-B-A/2 structure with thickness variation
M = M1*M2;
% M = M2*M1;
% M = M0*AA*M0; %%% defect structure (air)
Mm = M1*M2*M3;
% M = M2*M1*M4;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Reflection and Transmission
% R1(loop)=((ns*M(1,1)+ns*n0*M(1,2)-M(2,1)-n0*M(2,2))/(ns*M(1,1)+ns*n0*M(1,2)+M(2,1)+n0*M(2,2)));
% R(loop)=(R1(loop))*conj(R1(loop));
T1(loop)=(2*ns/(ns*M(1,1)+ns*n0*M(1,2)+M(2,1)+n0*M(2,2)));
T(loop)=(n0 / ns) * abs(T1(loop))^2;
Tt1(loop)=(2*ns/(ns*Mm(1,1)+ns*n0*Mm(1,2)+Mm(2,1)+n0*Mm(2,2)));
Tt(loop)=(n0 / ns) * abs(Tt1(loop))^2;
end
desired_freq = 60.895; % Change this to your desired peak frequency (THz)
peak_idx = find(f == desired_freq);
peak_val = Tt(peak_idx);
half_max = peak_val / 2;
peak_left_edge = find(Tt(1:peak_idx) < half_max, 1, 'last');
peak_right_edge = find(Tt(peak_idx:end) < half_max, 1, 'first') + peak_idx - 1;
peak_width = f(peak_right_edge) - f(peak_left_edge); %%% full width at half maximum
disp("peak width=")
peak width=
disp(peak_width)
0.2120
figure(1)
plot(f,T(:),'k') %%% Topological structure
hold on
plot(0.018+f,Tt(:),'r') %%% Topological coupled structure
hold off
xlim([60.5,61.5])
% xlim([58,67])
ylim([0,1])
xlabel("Frequency (THz)")
ylabel("Transmission")

请先登录，再回答此问题。

 采纳的回答

Mathieu NOE
Mathieu NOE 2024-3-22

0 个投票

hello
I simply added a for loop for the del parameter to see how peak amplitude and peak width of the red curve would evolve vs del (smooth or non smooth)
before that , I made a few mods to speed up the code
  • reduced the freq range to what you finally display (instead of computing lot more results that are then throw away)
  • modified how you do the peak width computation (here my suggestion will do linear interpolation, so we can use a slightly coarser freq vector without loss of accuracy)
if we use the new code with the complete range for del (from -1 to +1) , we got his result :
It's not evident to find a del value that correspond to a minimum point while having the transmission factor as close as possible to one. In fact it seems there are many possible values that correspond to this request
after some time , we can see that there is an interesting range between 0 and 0.2 that we can further investigate
narrowing down to the del 0 to 0.2 range gives this result
what would be the best del value here ?
seems to me del between 0.04 and 0.08 is a good value
now if you want to selct ONE optimum del value you have to define a criteria that combine peak amplitude and peak width with the correct weighting and minimize this value
I don't know right now how you would define the ponderation here - it's up to you
code :
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This is a program for 1-D Photonic crystal
% For Photonic band structure calculation and Transmission Spectrum
% This program uses Transfer Matrix Method
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc
clear
close all
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define material properties and simulation parameters
P = 0; %%% Hydrostatic pressure in MPa
n0 = 1; % Refractive index of air
ns = 1.46; % Refractive index of subtrate
n01 = 1.578; % Refractive index of PS
n02 = 1.484; % Refractive index of PMMA
% Calculate Epsilon for n01, n02
e01 = n01^2;
e02 = n02^2;
% Constants
p11_1 = 0.32;
p12_1 = 0.31;
p11_2 = 0.3;
p12_2 = 0.297;
v1 = 0.35;
E1 = 3.3e3; % Convert to scientific notation (3.3 x 10^3)
v2 = 0.37;
E2 = 3.0303e3; % Convert to scientific notation (3.0303 x 10^3)
% Formula for calculating nA
e1 = e01 - ((e01^2)/2) * (((p11_1/E1) * (v1+1) * P) + ((p12_1/E1) * (3*v1+1) * P));
nA = sqrt(e1);
% Formula for calculating nA
e2 = e02 - ((e02^2)/2) * (((p11_2/E2) * (v2+1) * P) + ((p12_2/E2) * (3*v2+1) * P));
nB = sqrt(e2);
% del_range= -1:0.01:1; %%% Thickness variation parameter (full range search)
del_range= 0:0.001:0.2; %%% Thickness variation parameter (narrow range search)
% Initialization of parameters
dair=1200e-9; % Thickness of air in meters
dA=780e-9; % Thickness of First Layer in meters
dB=830e-9; % Thickness of Second Layer in meters
% ddA=(1+del)*dA; % Thickness of First Layer with thickness variation in meters
% ddB=(1+del)*dB; % Thickness of Second Layer with thickness variation in meters
c=3e8; % Velocity of Light (m/s)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Frequency Range
% f=linspace(58,67,16001); % Frequncy in THz
% f=58:0.001:67; % Frequncy in THz
f=60.5:0.01:61.5; % Frequncy in THz
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for k = 1:numel(del_range)
del = del_range(k);
ddA=(1+del)*dA; % Thickness of First Layer with thickness variation in meters
ddB=(1+del)*dB; % Thickness of Second Layer with thickness variation in meters
for loop=1:numel(f)
%%% angular frequency
w=2*pi*f*1e12;
a=dA+dB; %%% Lattice Constant
DA=((w(loop))/c)*dA*nA; %%% wave number of first layer
DB=((w(loop))/c)*dB*nB; %%% wave number of second layer
Da=((w(loop))/c)*(dA*0.5)*nA; %%% wave number of first half layer
Db=((w(loop))/c)*(dB*0.5)*nB; %%% wave number of second half layer
DDA=((w(loop))/c)*ddA*nA; %%% wave number of first layer with thickness variation
DDB=((w(loop))/c)*ddB*nB; %%% wave number of second layer with thickness variation
DDa=((w(loop))/c)*(ddA*0.5)*nA; %%% wave number of first half layer with thickness variation
DDb=((w(loop))/c)*(ddB*0.5)*nB; %%% wave number of second half layer with thickness variation
Dair=((w(loop))/c)*(dair*0.5)*n0; %%% wave number of air layer
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Transfer Matrix elements of air
a11=cos(Dair); a12=-1i*sin(Dair)/n0; a21=-1i*n0*sin(Dair); a22=cos(Dair);
AA=[a11 a12; a21 a22];
%%% Transfer Matrix elements of first layer
m11=cos(DA); m12=-1i*sin(DA)/nA; m21=-1i*nA*sin(DA); m22=cos(DA);
mA=[m11 m12; m21 m22];
%%% Transfer MAtrix elements of Second layer
l11=cos(DB); l12=-1i*sin(DB)/nB; l21=-1i*nB*sin(DB); l22=cos(DB);
mB=[l11 l12; l21 l22];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Transfer Matrix elements of first half layer
ma11=cos(Da); ma12=-1i*sin(Da)/nA; ma21=-1i*nA*sin(Da); ma22=cos(Da);
ma=[ma11 ma12; ma21 ma22];
%%% Transfer MAtrix elements of Second half layer
lb11=cos(Db); lb12=-1i*sin(Db)/nB; lb21=-1i*nB*sin(Db); lb22=cos(Db);
mb=[lb11 lb12; lb21 lb22];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Transfer Matrix elements of first layer with thickness variation
dm11=cos(DDA); dm12=-1i*sin(DDA)/nA; dm21=-1i*nA*sin(DDA); dm22=cos(DDA);
mA1=[dm11 dm12; dm21 dm22];
%%% Transfer MAtrix elements of Second layer with thickness variation
dl11=cos(DDB); dl12=-1i*sin(DDB)/nB; dl21=-1i*nB*sin(DDB); dl22=cos(DDB);
mB2=[dl11 dl12; dl21 dl22];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Transfer Matrix elements of first half layer with thickness variation
dma11=cos(DDa); dma12=-1i*sin(DDa)/nA; dma21=-1i*nA*sin(DDa); dma22=cos(DDa);
ma1=[dma11 dma12; dma21 dma22];
%%% Transfer MAtrix elements of Second half layer with thickness variation
dlb11=cos(DDb); dlb12=-1i*sin(DDb)/nB; dlb21=-1i*nB*sin(DDb); dlb22=cos(DDb);
mb2=[dlb11 dlb12; dlb21 dlb22];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Unit cell
m=(mA*mB);
M0=(m^25); %%% A-B structure
M1=(mb*mA*mb)^25; %%% B/2-A-B/2 structure
M2=(ma*mB*ma)^25; %%% A/2-B-A/2 structure
M3=(mb2*mA1*mb2)^25; %%% B/2-A-B/2 structure with thickness variation
M4=(ma1*mB2*ma1)^25; %%% A/2-B-A/2 structure with thickness variation
M = M1*M2;
% M = M2*M1;
% M = M0*AA*M0; %%% defect structure (air)
Mm = M1*M2*M3;
% M = M2*M1*M4;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Reflection and Transmission
% R1(loop)=((ns*M(1,1)+ns*n0*M(1,2)-M(2,1)-n0*M(2,2))/(ns*M(1,1)+ns*n0*M(1,2)+M(2,1)+n0*M(2,2)));
% R(loop)=(R1(loop))*conj(R1(loop));
T1(loop)=(2*ns/(ns*M(1,1)+ns*n0*M(1,2)+M(2,1)+n0*M(2,2)));
T(loop)=(n0 / ns) * abs(T1(loop))^2;
Tt1(loop)=(2*ns/(ns*Mm(1,1)+ns*n0*Mm(1,2)+Mm(2,1)+n0*Mm(2,2)));
Tt(loop)=(n0 / ns) * abs(Tt1(loop))^2;
end
% red line data processing
ft = 0.018+f;
[peak_val(k),~] = max(Tt);
half_max = peak_val(k) / 2;
[f_left_edge,f_right_edge] = find_zc(f,T,half_max);
peak_width(k) = f_right_edge - f_left_edge;
end
plot(del_range,peak_val,'b',del_range,peak_width,'r');
xlabel("Del")
legend('peak val','peak width')
return
figure(1)
plot(f,T(:),'k') %%% Topological structure
hold on
plot(ft,Tt(:),'r') %%% Topological coupled structure
hold off
% xlim([60.5,61.5])
% xlim([58,67])
ylim([0,1])
xlabel("Frequency (THz)")
ylabel("Transmission")
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [ZxP,ZxN] = find_zc(x,y,threshold)
% put data in rows
x = x(:);
y = y(:);
% positive slope "zero" crossing detection, using linear interpolation
y = y - threshold;
zci = @(data) find(diff(sign(data))>0); %define function: returns indices of +ZCs
ix=zci(y); %find indices of + zero crossings of x
ZeroX = @(x0,y0,x1,y1) x0 - (y0.*(x0 - x1))./(y0 - y1); % Interpolated x value for Zero-Crossing
ZxP = ZeroX(x(ix),y(ix),x(ix+1),y(ix+1));
% negative slope "zero" crossing detection, using linear interpolation
zci = @(data) find(diff(sign(data))<0); %define function: returns indices of +ZCs
ix=zci(y); %find indices of + zero crossings of x
ZeroX = @(x0,y0,x1,y1) x0 - (y0.*(x0 - x1))./(y0 - y1); % Interpolated x value for Zero-Crossing
ZxN = ZeroX(x(ix),y(ix),x(ix+1),y(ix+1));
end

5 个评论
显示 3更早的评论 隐藏 3更早的评论

GULZAR
GULZAR 2024-3-26
I need a graph of transmission versus frequency corresponds to the optimized del_range
GULZAR
GULZAR 2024-3-26
Can you let me know
Which optimization technique this?
Mathieu NOE
Mathieu NOE 2024-3-26
to you first comment , it seems the last figure of my answer is almost what you need
simply pick the del value inside the 0 to 0.2 range that seems to best fit both of your goals (should be easy)
to your second comment : nothing fancy, simply created the "grid" - here a 1 D vector of del values , then computed your outputs vs this del range and find out which optimal values are best
some would refer to a "brute force approach" when you define all possible parameters values, compute the outputs and pick the best result. Of course , here we deal with a 1 D problem which moderate number of combinations, so it's quite fast
would be a different story with hundreds of parameters (dimensions)

请先登录，再进行评论。

更多回答（0 个）

类别

在 帮助中心 和 File Exchange 中查找有关 Oceanography and Hydrology 的更多信息

产品

版本

R2022a

标签

Community Treasure Hunt

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

Start Hunting!

Translated by