Grating coupler - Matlab-driven optimization (2D)代码解读
魏楷
Grating coupler - Matlab-driven optimization (2D)
Main函数
%添加Lumerical Matlab API
path(path,'C:\Program Files\Lumerical\FDTD\api\matlab'); %将Lumerical API路径添到Matlab；
sim_file_path=('C:\Learning\Matlab Nonlinear optimization'); %用matlab打开的fdtd文件所在的存放位置
sim_file_name=('grating_coupler_2D_Matlab_Optimization.fsp'); %用matlab打开的名为grating_coupler_2D_Matlab_Optimization.fsp的文件
%Open FDTD session
h=appopen('fdtd');
%Pass the path variables to FDTD
appputvar(h,'sim_file_path',sim_file_path); %通过h将Matlab工作空间中的变量sim_file_path作为变量sim_file_path发送到Lumerical的工作空间中
appputvar(h,'sim_file_name',sim_file_name);
%Load the FDTD simulation file and get simulation parameters
code=strcat('cd(sim_file_path);',...
'load(sim_file_name);',...
'select("grating_coupler_2D");',...
'coupler_y_pos=get("y");',... %strcat 水平串联字符串
'coupler_thickness=get("h total");',...
'select("fiber::source");',...
'source_y_pos=get("y");',...
'select("fiber");',...
'fiber_y_pos=get("y");',...
'gap=abs(fiber_y_pos+source_y_pos-coupler_y_pos-coupler_thickness);',...
'select("FDTD");',...
'FDTD_span=get("x span");');
%send the script in 'code' to Lumerical FDTD Solutions
appevalscript(h,code); %在fdtd中执行code命令
%Get variables 'FDTD_span' and 'gap' from FDTD workspace to Matlab
%workspace
global gap %定义全局变量gap
gap=appgetvar(h,'gap'); % 通过h从Lumerical工作空间中检索变量gap，并将其作为变量gap添加到Matlab工作空间中
FDTD_span=appgetvar(h,'FDTD_span');
%Function to be optimized x(x_position,theta,pitch,dc)
f=@(x,y)Coupler_Optimization(x(1),x(2),x(3),x(4),h); % @是定义句柄的运算符，相当于建立一个Coupler_Optimization函数文件，（x,y）表示函数文件里面的参数，多个参数用逗号隔开
%Optimization starting points, constraints and boundaries
%fmincon求解器中的输入参数
%The format is [x_position in um/10, theta in degrees/10, pitch in um, duty cycle]
x0=[0.5,1.7,0.75,0.75]; %起始点，指定为实向量或实数组。 求解器使用x0中的元素数量和大小来确定f接受的变量的数量和大小。sqp算法将超出边界的x0分量重置为相应边界的值。
%%对于大问题，可以将A、b、Aeq、beq作为稀疏向量传递
A=[]; %线性不等式约束，指定为实矩阵。 A是一个M×N的矩阵，其中M是不等式的数量，N是变量的数量（x0中 的元素数量）。
b=[]; %线性不等式约束，指定为实向量。 b是与A矩阵有关的M元素向量。 如果将b作为行向量传递，则求解器会在内部将b转换为列向量b。
Aeq = []; %线性相等约束，指定为实矩阵。 Aeq是一个Me-by-N矩阵，其中Me是等式数，N是变量数（x0中的元素数）。
beq=[]; %线性相等约束，指定为实向量。 beq是与Aeq矩阵有关的Me元素向量。 如果将beq作为行向量传递，则求解器会在内部将beq转换为列向量beq。
lb=[0,1.1,0.5,0.5]; %下限，所有的x(i)>=lb(i)，lb中向量个数必须少于x0中向量个数，否则发出警告
ub=[FDTD_span/2e-5,2,0.8,0.8]; %上限，所有的x(i)<=ub(i)，ub中向量个数必须少于x0中向量个数，否则发出警告
nonlcon=@confun; %建立一个confun函数文件
nonlcon
nonlcon是一个接受向量x或数组x并返回两个数组c（x）和ceq（x）的函数。
c（x）是x处的非线性不等式约束的数组。 fmincon试图满足对于c的所有条目，c（x）<= 0。
ceq（x）是x处的非线性等式约束的数组。 fmincon试图满足对于ceq的所有条目，ceq（x）= 0。
%Optimization settings
%fmincon多变量非线性约束优化算法，发现满足约束的局部最小值
options = optimoptions('fmincon'); %返回fmincon求解器的一组默认选项
%设置fmincon求解器的属性
options = optimoptions(options,'FiniteDifferenceStepSize',1e-3); %标量或矢量步长因子有限差分
options = optimoptions(options,'Algorithm','sqp'); %alt: interior-point，选择sqp优化算法
options = optimoptions(options,'MaxIter', 100); %允许的最大迭代次数
options = optimoptions(options,'PlotFcns', { @optimplotfval }); %在算法执行时绘制各种进度度量，对于自定义绘图函数，传递函数句柄
options = optimoptions(options,'Display','iter-detailed'); %显示级别
options = optimoptions(options,'StepTolerance',1e-3); %步进公差
%Run optimization
[x,fval]=fmincon(f,x0,A,b,Aeq,beq,lb,ub,nonlcon,options); %返回解x处目标函数f的值
T_avg=1-fval;
FiniteDifferenceStepSize
当FiniteDifferenceStepSize设置为向量v时，正向有限差分量为
delta = v.*sign′(x).*max(abs(x),TypicalX);
其中sign′(x) = sign(x)， sign′(0) = 1除外。
中心有限差分为
delta = v.*max(abs(x),TypicalX);
标量FiniteDifferenceStepSize扩展为向量时， 对于正向有限差分，默认值为sqrt（eps），对于中央有限差分，默认值为eps ^（1/3）。
fmincon Algorithms
-
‘interior-point’(default)
处理大型的，稀疏的问题以及小型的密集问题，在所有迭代中都满足边界，是一种大规模算法。 -
‘trust-region-reflective’
要求提供一个渐变，并且仅允许边界或线性相等约束，而不能同时提供两者。在这些限制内，该算法可以有效地处理大型稀疏问题和小型密集问题。这是一种大规模算法。 -
‘sqp’
在所有迭代中都满足边界，不是大规模算法。 -
‘sqp-legacy’
与’sqp’类似，但是通常更慢并且使用更多的内存。 -
‘active-set’
可以采取较大的步骤，从而提高速度。该算法在约束不平滑的某些问题上有效。这不是大规模算法。
PlotFcn
‘optimplotx’ 绘制当前点
‘optimplotfunccount’ 绘制函数计数
‘optimplotfval’ 绘制函数值
‘optimplotfvalconstr’ 将找到的最佳可行目标函数值绘制为折线图。该图使用1e-6的可行性公差将不可行点显示为红色，将可行点显示为蓝色
‘optimplotconstrviolation’ 绘制最大约束违例
‘optimplotstepsize’ 绘制步长
‘optimplotfirstorderopt’ 绘制一阶最优度量
Display
‘off’ or ‘none’ 不显示任何输出。
‘iter’ 在每次迭代时显示输出，并给出默认的退出消息。
‘iter-detailed’ 在每次迭代时显示输出，并给出技术退出消息。
‘notify’ 仅在函数未收敛时显示输出，并提供默认退出消息。
‘notify-detailed’ 仅在功能未收敛时显示输出，并给出技术退出消息。
‘final’ (default) 仅显示最终输出，并提供默认退出消息。
‘final-detailed’ 仅显示最终输出，并给出技术退出消息。
%Display optimization results
disp(strcat({'Optimized x position: '},num2str(x(1)*10),{' um'})); %strcat水平串联字符串
disp(strcat({'Optimized angle theta: '},num2str(x(2)*10),{' degrees'})); %num2str将括号里面的内容转换为字符串
disp(strcat({'Optimized pitch: '},num2str(x(3)),{' um'})); %disp %disp显示变量的值
disp(strcat({'Optimized Duty cycle: '},num2str(x(4))));
disp(strcat({'Optimized T_avg: '},num2str(T_avg)));
%Close session
appclose(h);
Coupler_Optimization
function y=Coupler_Optimization(x_position,theta,pitch,dc,h) %需要优化的4个参数，即输入参数
%Modify fiber position, theta, duty cycle, pitch and run the
%simulation
code=strcat('switchtolayout;',...
'select("grating_coupler_2D");',...
'set("pitch",',num2str(pitch*1e-6,16),');',... %num2str(x,精度)将x转换为逗号后面对应精度形式的数
'set("duty cycle",',num2str(dc,16),');',...
'select("fiber");',...
'set("x",',num2str(x_position*1e-5,16),');',...
'set("theta0",',num2str(theta*10,16),');',...
'run;');
appevalscript(h,code); %在fdtd中执行code命令
%Get the coupled power from T monitor to
%FDTD workspace as variable 'T_avg_FDTD'
code=strcat('T_avg_FDTD=getresult("model","T_avg");');
appevalscript(h,code); %fdtd仿真运行后得到T监视器的T_avg值，赋给code
%Get the average transmission(figure of merit) from FDTD workspace to
%Matlab workspace
T_MatlabFun=appgetvar(h,'T_avg_FDTD'); % 通过h从Lumerical工作空间中检索变量T_avg_FDTD，并将其作为变量T_MatlabFun添加到Matlab工作空间中
y=1-abs(T_MatlabFun); %返回1-|T_MatlabFun|
%Uncommnet this section to display the optimized parameters and
%figure of merit for each run in FDTD
%disp(strcat('x_pos= ',num2str(x_position,10),'...theta= ',num2str(theta,10),'...pitch= ',num2str(pitch,10),'...DC= ',num2str(dc,10)));
%disp(strcat('1-abs(T)= ',num2str(y,10)));
end
confun
function [c, ceq] = confun(x)
% Nonlinear inequality constraints x(x_pos,theta,pitch,dc)
global gap
c = tan(x(2)*pi/180)*gap-x(1)*1e-5; %
% Nonlinear equality constraints
ceq = [];
%%优化的目的是在1500nm至1600nm的波长范围内最大化进入SOI波导模式的平均传输
类似资料：