Filtered gradient with external optimizer (IPOPT)


I would like to perform some topology optimization using the IPOPT optimizer and the automatic differentiation functionalities provided by dolfin-adjoint.

However, I would like the optimizer to account for a filtered/modified version of the gradient of my functional and not directly the one resulting from the operations saved on the Tape. This is referred to as “sensitivity-filtering” in the topology optimization field.

Do you have suggestions or advice on how to pass my custom gradient to the solver? And how can it take this filtering operation into account for every iteration of the optimization?

Here is a MWE based on the tutorial for topology optimization of heat conduction problems. I have indicated the default gradient and the desired one.

from __future__ import print_function
from dolfin import *
from dolfin_adjoint import *

from pyadjoint.drivers import compute_gradient
import matplotlib.pyplot as plt

import numpy as np
#import sklearn.metrics.pairwise as sp

    from pyadjoint import ipopt
except ImportError:
    info_red("""This example depends on IPOPT and pyipopt. \
  When compiling IPOPT, make sure to link against HSL, as it \
  is a necessity for practical problems.""")

parameters["std_out_all_processes"] = True

p = Constant(5)        # power used in the solid isotropic material with penalisation (SIMP) rule
eps = Constant(1.0e-3) # epsilon used in the solid isotropic material
rmin = 0.05 # filter radius
alpha = Constant(1.0e-8) # regularisation coefficient in functional
V_star = Constant(0.4)      # volume bound on the control

def k(a):
    """Solid isotropic material with penalisation (SIMP) conductivity
  rule, equation (11)."""
    return eps + (1 - eps) * a**p

n = 100
mesh = UnitSquareMesh(n, n)
A = FunctionSpace(mesh, "CG", 1)  # function space for control
P = FunctionSpace(mesh, "CG", 1)  # function space for solution

class WestNorth(SubDomain):
    """The top and left boundary of the unitsquare, used to enforce the Dirichlet boundary condition."""
    def inside(self, x, on_boundary):
        return (x[0] == 0.0 or x[1] == 1.0) and on_boundary

bc = [DirichletBC(P, 0.0, WestNorth())]
f = interpolate(Constant(1.0e-2), P) # the volume source term for the PDE

def forward(a):
    """Solve the forward problem for a given material distribution a(x)."""
    T = Function(P, name="Temperature")
    v = TestFunction(P)
    F = inner(grad(v), k(a)*grad(T))*dx - f*v*dx
    solve(F == 0, T, bc, solver_parameters={"newton_solver": {"absolute_tolerance": 1.0e-7, "maximum_iterations": 20}})
    return T

def J(a,T):
    """Functional of the optimization problem, with `a` control var and `T` temperature."""
    return assemble(f*T*dx + alpha * inner(grad(a), grad(a))*dx)

def filter_sensitivity(s_ini,radius_f,FS):
    """PDE Helmholtz-type filter for sensitivity"""
    v_f = TrialFunction(FS)
    w_f = TestFunction(FS)
    sys_left = (radius_f**2)*inner(grad(v_f), grad(w_f))*dx +v_f*w_f*dx
    sys_right = s_ini*w_f*dx
    sys_mat, sys_vec = assemble_system(sys_left,sys_right)
    s_f = Function(FS)
    solve(sys_mat, s_f.vector(), sys_vec, annotate=True)
    return s_f

if __name__ == "__main__":

    a = interpolate(V_star, A)
    T = forward(a)        

    Jfunc = J(a,T)
    m = Control(a)
    Jhat = ReducedFunctional(Jfunc, m)
    dJda = compute_gradient(Jfunc,Control(a)) # default gradient used to solve the minimization problem "behind the scenes"
    dJda_f = filter_sensitivity(dJda,rmin,a.function_space()) # desired gradient (to be performed every iteration)
    volume_constraint = UFLInequalityConstraint((V_star - a)*dx, m)

    problem = MinimizationProblem(Jhat, bounds=(0.0, 1.0), constraints=volume_constraint)

    solver = IPOPTSolver(problem, parameters={"acceptable_tol": 1.0e-3, "maximum_iterations": 100})
    a_opt = solver.solve()

    # Visualization
    Jfunc_opt = J(a_opt,forward(a_opt))
    dJda_opt = compute_gradient(Jfunc_opt,Control(a_opt))
    dJda_opt_f = filter_sensitivity(dJda_opt,rmin,a_opt.function_space())
    pl, ax = plt.subplots(); fig = plt.gcf(); fig.set_size_inches(12, 4)
    plt.subplot(1,3,1); p = plot(a_opt, title='Topology', mode='color');cbar = plt.colorbar(p)
    plt.subplot(1,3,2); p = plot(dJda_opt, title='Sensitivity', mode='color'); cbar = plt.colorbar(p) 
    plt.subplot(1,3,3); p = plot(dJda_opt_f, title='Sensitivity filtered', mode='color'); cbar = plt.colorbar(p) 

Thank you for your time.


An alternative to your problem would be to implement your own “optimizer” such as the Optimality Criteria method. This has been done in this paper (with the FEniCS code).

However, I don’t know either how to pass a custom gradient to existing solvers such as IPOPT …

Dear Brad,

Thank you for your answer and suggestion.
However, I would like to be able to use the IPOPT solver in this case.

Does anyone have more information regarding this?