# Filtered gradient with external optimizer (IPOPT)

Hello,

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 *

import matplotlib.pyplot as plt

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

try:
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.""")
raise

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)
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."""

"""PDE Helmholtz-type filter for sensitivity"""
v_f = TrialFunction(FS)
w_f = TestFunction(FS)

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_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)

``````

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 would like to be able to use the `IPOPT` solver in this case.