optimal_control.py

""" Solves the mother problem in optimal control of PDEs """

from dolfin import *
from dolfin_adjoint import *

dolfin.set_log_level(ERROR)
dolfin.parameters["optimization"]["test_gradient"] = False

# Define the domain [-1, 1] x [-1, 1]
n = 200
mesh = RectangleMesh(-1, -1, 1, 1, n, n)
# Define the function spaces
V = FunctionSpace(mesh, "CG", degree = 1)
W = FunctionSpace(mesh, "DG", degree = 0)
# Define the functions
u = Function(V, name = "Solution")
m = Function(W, name = "Control")
v = TestFunction(V)

def save_plot(x, name):
    filename = "_".join(name.split()).lower()
    try:
        xcg = project(x, V, annotate = False)
        viz_x = plot(xcg,
                wireframe = False,
                title = name,
                rescale = True,
                axes = False,
                basename = name,
                )
        #viz_x.elevate(-65) # tilt camera
        try:
            viz_x.set_min_max(0, 0.5*x.vector().max()) # color scale
        except:
            pass
        viz_x.update(xcg)
        viz_x.write_png(filename)
        viz_x.write_pdf(filename)

    except:
        pass
    
    out = File(filename+".pvd")
    out << project(x, W, annotate = False)

# Solve the forward model to create the tape 
F = (inner(grad(u), grad(v)) - m*v)*dx 
bc = DirichletBC(V, 0.0, "on_boundary")
solve(F == 0, u, bc)

x = triangle.x
# Define the functional of interest
u_d = exp(-1/(1-x[0]*x[0])-1/(1-x[1]*x[1]))
J = Functional((0.5*inner(u-u_d, u-u_d))*dx)
u_dp = project(u_d, V)
save_plot(u_dp, "Desired state")

# Define the reduced functional
J_hat = ReducedFunctional(J, SteadyParameter(m))
# Solve the optimisation problem
m_opt = minimize(J_hat, method = "L-BFGS-B", bounds = (0.0, 0.5), 
        options = {"gtol": 1e-16, "ftol": 1e-16, "maxiter": 10})

save_plot(m_opt, "Optimal control")
F_opt = replace(F, {m: m_opt})
solve(F_opt == 0, u, bc)
save_plot(u, "Optimal state")
save_plot(u-u_d, "Difference between optimal and desired state")