Non-density based adjoint optimization

[kalosu]

Hello there friends from the meep community,

I wanted to ask for a bit of your guidance. I have read the publication "High-performance hybrid time/frequency-domain topology optimization for large-scale photonics inverse design" in where you combine a density based design approach with adjoint problems based gradient evaluations.

I have no experience with such optimization methods and therefore I wanted to ask for your feedback about a similar but probably different implementation. Instead of an optimization based on individual pixels, I would like to optimize for a sufficiently "smooth" structure that could be described maybe by a linear combination of basis terms or an aperiodic grating structure for which each local grating period could be adjusted.

Since I imagine (again I am not really sure) that performing such an optimization would be more efficient through gradient-based methods, I would like to know if for such a case I could compute the gradients by relying on a similar adjoint problem formulation.

In summary, I would like to use a gradient-based optimization algorithm (via adjoint methods) but for "regular" and parametrizable structures with probably a large number of parameters (i.e, smooth surfaces and aperiod gratings)

Therefore, my questions would be:

1. Is what I am intending to do feasible?
2. If so, could I use the same code base from your examples?
   2.1 If not, what would I need to adjust from your code base?
3. In case the answer to 1 is no, what else would you recommend?

Thanks a lot for your feedback!

[smartalecH]

Check out some of the comments from #2067.

Basically what you are describing requires a level set (an "explicit" level set in this case). Luckily meep's adjoint code supports that (using the new subpixel smoothing feature), but you still need to map whatever basis functions you specify to densities. Again, the above discussion describes that.

We've been working on a paper that describes this for quite some time... life just gets in the way.

[kalosu]

Hi @smartalecH , thanks for the feedback and nice info.

So if I understand it correctly (after reading through the posts) I would need to find a mapping function that allows me to describe my structure as a function of the parameters over which I want to optimize. Such a mapping function should generate a required levelset used to compute material density distributions within MEEP. Is this understanding correct?

Also, I have read that there are some constraints that need to be satisfied:

1. The parameters of interest are "mapped" to an underlying grid (essentially a density-grid). For example, your "mapping function" could simply look like ====> What I tried to describe above

2. This function mapping must be "smooth" (you can't "threshold" or "project" yet). ====> So this means that I cannot have a structure with vertical sidewalls as in an aperiodic grating? Or does this mean that the mapping function itself should be smooth resulting in sharp material transitions that is mapped to a well defined value from the levelset?

3. the MaterialGrid's beta parameter is set to beta=np.inf. In other words, the projection happens "under the hood". ====> What do you mean here by the projection happens under the hood?

If you can provide me with some further feedback about these points it would be nice!

In general I would like to be able to perform an optimization for a structure that has well defined sharp edges. (Of course to fabricate structures with sharp and clear vertical edges is not always possible but I would prefer to then evaluate the performance of an "ideal" structure once I introduce such "imperfections" into the design).

Again thanks for your feedback!

[smartalecH]

Or does this mean that the mapping function itself should be smooth

This. By projecting to infinity, you'll create a discontinuous structure (like you're looking for). But in order for the subpixel smoothing to work, the function itself must be smooth (we use a first order approximation, and that breaks down if it's not smooth).

What do you mean here by the projection happens under the hood?

In the tutorials, we manually project inside the mapping function. Don't do that anymore. Rather, the MaterialGrid object (and also the OptimizationProblem.__call__() method) have a beta parameter that you use instead. The actual projection will happen in the C++ code (and do the averaging with it).

[kalosu]

Hi again @smartalecH, thanks for the clarification.

From what I understand, for a MaterialGrid object the pixel points defined within the grid can take only one of two possible material values which are specified as parameters to the MaterialGrid object.

Also, if I understand it correctly, the mapping function that I end up using should return a scalar value evaluated at each coordinate of my design region which depends on my optimization parameters and that should be bounded between 0 and 1, like in the example you have provided before:

import autograd.numpy as npa
x = npa.linspace(0,L,N) # where N is the number of MaterialGrid DOF, and L is grating length
def fx(Λ):
    return npa.sin(2*npa.pi/Λ*x) # where Λ is grating period

Nevertheless, each individual pixel from the MaterialGrid object will still be populated by one of the two material values right? (I guess this is what you mean by the internal projection with $\beta$ = $\infty$).

Based on this, I have the following questions:

1. Is my interpretation that the mapping function should return a value between 0 and 1 for each point in space correct?
2. How can I guarantee that the medium properties within a solid volume will not have discontinuities (let's say I just want to optimize the individual width of multiple rectangular blocks but still preserve the rectangular nature of each block)

Again, thanks a lot for your feedback!

[smartalecH]

Is my interpretation that the mapping function should return a value between 0 and 1 for each point in space correct?

Correct.

Edit: maybe not each point in space, but each point in your material grid.

How can I guarantee that the medium properties within a solid volume will not have discontinuities (let's say I just want to optimize the individual width of multiple rectangular blocks but still preserve the rectangular nature of each block)

This is where your mapping function is important. If you use a sin function like you have there, then after the projection you will always have rectangular slabs.

But I think what your asking is you want to have a "grid" of blocks, and you want to optimize their widths? I'm not sure I understand the question...

[stevengj]

In general I would like to be able to perform an optimization for a structure that has well defined sharp edges.

This happens anyway with density-based topology optimization once β goes to infinity (and one can also penalize intermediate values in other ways, e.g. with the damping parameter or a penalty term). (This kind of "binarization" process is usually necessary to obtain a manufacturable structure!) So it's not clear to me what new thing you're trying to accomplish here.

Yes, it's possible to implement level-set methods with custom degrees of freedom if you define a smooth mapping from your degrees of freedom to a level-set function with β=∞, as @smartalecH outlined above, but that's a lot of work to go through if you can use density-based TO as-is.

(Note that even with density-based TO you can still restrict the dimensionality of the design. e.g. if you want to design aperiodic gratings, i.e. univariate materials, you can just apply a 1d material grid to a 2d or 3d object. We do this all the time to design 2d-patterned structures in 3d, by applying a 2d material grid to a 3d design region.)

[kalosu]

Hi, thanks a lot for the feedback.

Actually, one of the reasons why I want to try running my optimization as a function of parameters different than the point-wise material properties is because in this form I could also have some insight behind the effect of these parameters into my design's performance.

I can see this from the paper "Inverse design of near unity efficiency perfectly vertical grating couplers" (https://opg.optica.org/directpdfaccess/3bbbc34a-0307-44c2-b9b11d92e5349e6b_381702/oe-26-4-4766.pdf?da=1&id=381702&seq=0&mobile=no) in where the local grating periods for a fixed number of grating "unit cells" are used as the parameters varied over the optimization. I feel that it is a bit more "intuitive" to evaluate the influence of each single parameter on the overall performance in this form.

In any case, do you think is it worth trying to use the adjoint method for such a case? Would it not be better to implement this using a gradient free optimization routine instead? (As far as I understood it, the adjoint method based optimization approach is better suited for large parameter spaces such as in the case of having individual material pixels as degrees of freedom)

[stevengj]

Gradient-based optimization (via adjoint-method gradients) is normally strictly faster than derivative-free optimization. However, if you have only a few parameters (probably < 10), derivative-free methods are feasible and have the advantage of being extremely easy to use for whatever parameters you wish.

[kalosu]

Hi @smartalecH , hi @stevengj ,

Thanks a lot for your help and feedback. I have given it a try first with a "simple" case.

I have optimized an aspheric surface (using the conic constant and the radius of curvature as parameters) via another software package.

For such a surface, I define a uniform illumination from air for which the aspheric surface focuses the incoming field. Following, I compute the coupling efficiency or the similarity of the focused field with respect to a Gaussian like field distribution.

I have implemented the same structure and system using meep for which I verified that I obtain similar results as for the reference design software.

I then decided to "adjust" the parameters describing such a surface in order to give it a try to the adjoint based optimization approach. (I know that I am optimizing just for two simple parameters but I wanted to start by optimizing a "simple" reference system with few parameters)

The problem that I have is that the merit function does not change so much in value and after a few iterations, the optimization routine just finishes without having the merit function reaching to the value from the original design. The following curve shows an example of the evolution of my merit function:

Also, I have noticed that the gradient with respect to my two free parameters is really "small" (in the order of 10e-6)

Could it be that I am not providing the correct gradient information to the optimization routine?

This is the mapping function that I am using for generating the MaterialGrid distribution:

def fx(rc: float, kappa: float):
	xcoord = npa.linspace(-design_region_width*0.5,design_region_width*0.5,Nx) 
	ycoord = npa.linspace(design_region_height_center-design_region_height*0.5,design_region_height_center+design_region_height*0.5,Ny)   
	xv, yv = npa.meshgrid(xcoord,ycoord)
	c = 1 / (rc) 

	zz = -((0 + ((c*npa.power(xv,2)))/(1+npa.sqrt(1-(1+kappa)*(npa.power(c,2))*(npa.power(xv,2)))))) 

	mapping = npa.where(zz-yv>=0,npa.where(npa.abs(xv)<=design_region_width*0.5,1/(1+npa.exp(-1*(-yv-zz))),0),0      ) 

	mapping = mapping.flatten(order='F') 

	return mapping

For the mapping, I am using a sigmoid function with the threshold value defined at the contour of the aspheric surface height profile (zz). As I said, I have verified that above mapping function results in the reference aspheric surface if I directly use the parameters obtained during the original optimization via the other commercial software.

As for the function that is given to the optimization routine, this is how I have defined it:

def f(opti_var,gradient):
	print ("Current iteration: {}".format(cur_iter[0]+1))
	print ("Current cr value: {}".format(opti_var[0]))
	print ("Current kappa: {}".format(opti_var[1]))

	f,dJ_du = opt([fx(opti_var[0],opti_var[1])])
	if gradient.size>0:
		gradient[0] = tensor_jacobian_product(fx,0)(opti_var[0],opti_var[1],dJ_du)
		gradient[1] = tensor_jacobian_product(fx,1)(opti_var[0],opti_var[1],dJ_du)

	evaluation_history.append(npa.real(f))
	cur_iter[0] = cur_iter[0] +1
	print (f"Objective value: {f}")
	print (f"dJ_du: {dJ_du}")
	print (f"gradients: {gradient}")
	return np.real(f)

As you can see, I use the ```tensor_jacobian_product" function in order to back propagate the gradients from the pointwise gradient evaluations for each pixel point to the gradients corresponding to my two optimization parameters.

Finally, this is how I set up everything in order to run the optimization routine:

algorithm = nlopt.LD_MMA 
lb = np.asarray([30.0,-1]) 
ub = np.asarray([60.0,1.0]) 

solver = nlopt.opt(algorithm,2) 
solver.set_lower_bounds(lb)
solver.set_upper_bounds(ub)
solver.set_max_objective(lambda a,g: f(a,g))
solver.set_maxeval(200)
solver.set_ftol_rel(1e-5)
xopt = solver.optimize(x_global) 
print ("Optimization finished!!")

I have used different MEEP example scripts as references and adapted my code accordingly as shown above. Do you notice anything suspicious or problematic in how I set up everything?

Also, do you have any idea in mind about how could I tackle this problem? What could I test? Could it be that my mapping function is not adequate? Should I try a different optimization algorithm or should I use a different formulation for the optimization problem? (Using an epigraph formulation?) Any suggestions on how could I move on forward with this?

As always, thanks a lot for your feedback and help!