In this example, we will demonstrate how MATLAB can be used to drive a multivariable nonlinear optimization of a grating coupler in FDTD via Lumerical's Automation API. In addition, we will demonstrate how to setup a MATLAB function based on arbitrary simulation parameters to specify a nonlinear constraint for the optimization.
While we focus on a specific example of a grating coupler in FDTD, the target audience includes anyone who is interested in driving and controlling Lumerical applications directly from MATLAB.
Requirements
MATLAB with Optimization Toolbox
MATLAB API Setup and Workflow
The key component required to drive simulations from MATLAB is the automation API that comes with 2016B or newer versions of Lumerical products. To allow for correct functionality, OS must know the location of Lumerical applications. On Windows machines this setup is executed automatically during the installation process, but it might be necessary for Linux users to add the install directory of given Lumerical product to their system path as described on this page.
The second step required for smooth functionality involves telling MATLAB where to look for the API before using any interoperability related commands. The MATLAB API is located in the Lumerical installation folder. In our example when we use FDTD on a Windows machine this would be typically under the following location:
C:\Program Files\Lumerical\FDTD\api\matlab
One option is to run the following command every time after MATLAB is launched and prior to using API related commands:
path(path,'C:\Program Files\Lumerical\FDTD\api\matlab');
The other, more convenient option is to add the path permanently to the startup.m file that is loaded every time MATLAB is launched.
Once the MATLAB API setup is completed we can use the related interoperability commands referenced in the "See Also" section. The basic API workflow is depicted in figure 1 and it consists of four operations:
 Open one or more Lumerical session
 Send and execute Lumerical script commands to one of the active sessions
 Transfer variables between the Lumerical and the MATLAB workspace via an active session
 Close Lumerical session
Note that operations 2,3 and 4 can be executed only by referencing to an active session otherwise they will result in an error. The session referencing is done by using an unique ID that is assigned to each session during its initiation. This ID system allows the users to open and control multiple sessions at the same time.
Optimization Project Goal
For this optimization project we use the FDTD simulation file from the Grating Coupler 2DFDTD example referenced above. This will allow us to compare the results from the MATLAB optimization with the results obtained by using a combination of Lumerical's builtin parameter sweep and particle swarm optimization utility. The goal of the optimization is to maximize the average transmission into the SOI waveguide mode in the wavelength range of 1500nm to 1600nm.
For this example, the following four parameters are chosen as the optimization variables:
 Angle theta representing the source injection angle in respect to y axis
 Source position in respect to the grating in x direction
 Grating pitch
 Grating duty cycle
In addition, this example demonstrates how to use a nonlinear constraint function that requires a specific FDTD simulation parameter as an input in addition to the variables used for optimization. To do this, we will constrain the x position of the source (x_{s}) so the beam axis at the surface of the grating does not extend beyond the end of the grating itself (x_{0}). This cannot be handled by a simple boundary condition since both angle theta and x position of the source are changing during the optimization. Moreover, we will use the MATLAB API to obtain the information about the distance between the source and the surface of the grating, both being inputs to the constraint function as follows:
$$x_s \geq \mathrm{tan}(\theta)\times gap$$
Figure 3: Visualization of the nonlinear constraint related to the source position and angle theta.
Optimization Project Setup
This project takes advantage of MATLAB "fmincon" multivariable nonlinear constrained optimization algorithm and it consists of four files:
grating_coupler_2D_MATLAB_Optimization.fsp 
2D FDTD simulation file. This file is launched from the main MATLAB script via automation API and it returns the average transmission as a function of the input parameters provided by the optimization script 
Coupler_Optimization_Main_Script.m 
Main MATLAB script file used to:

Coupler_Optimization.m 
MATLAB function that is called during the optimization. It receives set of parameters from fmincon and sets up the FDTD simulation accordingly. Once the simulation is completed in FDTD, this function extracts the average transmission and passes it back to the optimization routine. Note that since fmincon searches for the global minimum, we specify the figure of merit as fval=1T_avg; 
confun.m 
Nonlinear constraint function. Its variables are angle theta, x position of the source and the gap between the source and the grating as shown on figure 3. 
To run the example, download all three MATLAB files and the FDTD simulation file into the same folder. Open Coupler_Optimization_Main_Script.m and modify the sim_file_path to the path for the folder containing all downloaded simulation files to specify the location of the FDTD simulation file. Once the path is updated, simply run the simulation file and wait until MATLAB runs multiple optimization iterations. The final values of the optimized parameters and average transmission will be displayed in the MATLAB command window. The MATLAB script files are commented for clarity and can be reviewed in the section below:
Note: Prior running the example, you need to update the main script for sim_file_path to reflect your folder structure 
Coupler_Optimization_Main_Script.m
clear; %Add Lumerical Matlab API path path(path,'C:\Program Files\Lumerical\2019b\api\matlab'); sim_file_path=('Change\This\To\Your\Folder\Structure'); % update this path to user's folder sim_file_name=('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); 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");',... '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_poscoupler_y_poscoupler_thickness);',... 'select("FDTD");',... 'FDTD_span=get("x span");'); %send the script in 'code' to Lumerical FDTD appevalscript(h,code); %Get variables 'FDTD_span' and 'gap' from FDTD workspace to Matlab %workspace global gap gap=appgetvar(h,'gap'); 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); %Optimization starting points, constraints and boundaries %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]; A=[]; b=[]; Aeq = []; beq=[]; lb=[0,1.1,0.5,0.5]; ub=[FDTD_span/2e5,2,0.8,0.8]; nonlcon=@confun; %Optimization settings options = optimoptions('fmincon'); options = optimoptions(options,'FiniteDifferenceStepSize',1e3); options = optimoptions(options,'Algorithm','sqp'); %alt: interiorpoint options = optimoptions(options,'MaxIter', 100); options = optimoptions(options,'PlotFcns', { @optimplotfval }); options = optimoptions(options,'Display','iterdetailed'); options = optimoptions(options,'StepTolerance',1e3); %Run optimization [x,fval]=fmincon(f,x0,A,b,Aeq,beq,lb,ub,nonlcon,options); T_avg=1fval; %Display optimization results disp(strcat({'Optimized x position: '},num2str(x(1)*10),{' um'})); disp(strcat({'Optimized angle theta: '},num2str(x(2)*10),{' degrees'})); disp(strcat({'Optimized pitch: '},num2str(x(3)),{' um'})); disp(strcat({'Optimized Duty cycle: '},num2str(x(4)))); disp(strcat({'Optimized T_avg: '},num2str(T_avg))); %Close session appclose(h);
Coupler_Optimization.m
function y=Coupler_Optimization(x_position,theta,pitch,dc,h) %Modify fiber position, theta, duty cycle, pitch and run the %simulation code=strcat('switchtolayout;',... 'select("grating_coupler_2D");',... 'set("pitch",',num2str(pitch*1e6,16),');',... 'set("duty cycle",',num2str(dc,16),');',... 'select("fiber");',... 'set("x",',num2str(x_position*1e5,16),');',... 'set("theta0",',num2str(theta*10,16),');',... 'run;'); appevalscript(h,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); %Get the average transmission(figure of merit) from FDTD workspace to %MATLAB workspace T_MATLABFun=appgetvar(h,'T_avg_FDTD'); y=1abs(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('1abs(T)= ',num2str(y,10))); end
confun.m
function [c, ceq] = confun(x) % Nonlinear inequality constraints x(x_pos,theta,pitch,dc) global gap c = tan(x(2)*pi/180)*gapx(1)*1e5; % Nonlinear equality constraints is left empty ceq = [];
Optimization Results
The maximum average transmission achieved with the MATLAB driven optimization is ~40%, which is in good agreement with the value obtained using the Lumerical builtin parameter sweep/particle swarm optimization routines. The optimized values for all parameters shown in table 1 are close to the reference example demonstrating that we are able to reliably find the optimal parameters using this methodology. Note that the exact values might slightly change in each optimization run as the function minimum can be achieved with a range of parameter values rather than with one exact set of values.
Parameter 
MATLAB optimization value 
Reference example value 
Theta [degrees] 
12.21 
13.45 (fixed) 
Source x position [um] 
3.15 
3 um (from parameter sweep) 
Pitch [um] 
0.6475 
0.66 
Duty Cycle 
0.67 
0.66 
Optimized average transmission 
0.394 
0.388 
Table 1: Comparison between the optimized parameters and the reference example
The time required for the fourvariable optimization is comparable to the time required for the twovariable optimization with the particle swarm method, and it can be considerably influenced by setting the starting points, boundaries and the optimization parameters. Figure 4 shows comparison of two identical optimization setups that differ only in the finite difference step size (FDSS). It is clear that smaller step takes longer to estimate the trends and close in on the optimal values. As a result, we were able to achieve the optimized transmission within 8 iterations with FDSS=1e3, while 14 iterations were required with FDSS=1e4.
Figure 4: Tracking of the figure of merit(fval=1Average Transmission) for different finite difference step size (left shows 1e3 and right shows 1e4)
Optimization tips
 If the finite difference step size is too small, the simulation result might not change significantly and the optimization algorithm might terminate prematurely. For example, in our 2D FDTD simulation, moving the source by 1nm in x direction will have minimal to no impact on the resulting average transmission.
 It is a good idea to run the optimization multiple times with slightly different starting points and boundaries to verify that the results are comparable.
 If the optimization ends with one or more parameters having the boundary value, it is likely that the minimum is artificial and the boundary should be moved.
 Sometimes it can be useful to have the optimization parameters within the same order of magnitude.