optimal_control_transient.py

''' Example for solving a nonlinear optimisation problem with time-dependent control.

It solves:

min_{u,m} 1/2 ||u(t=T)-u_d(t=T)||^2 

subject to:

\grad u/\grad t - \nabla^2 u + u^3 = m,

where m: \Omega X [0, T] -> R

'''
from dolfin import *
from dolfin_adjoint import *
dolfin.set_log_level(ERROR)

# Define the discrete function(spaces)
N = 4 # Number of timesteps
timestep = 0.01
n = 50 # Mesh refinement
mesh = RectangleMesh(-1, -1, 1, 1, n, n)
V = FunctionSpace(mesh, "CG", 1)
W = FunctionSpace(mesh, "DG", 0)
u = Function(V, name="State")
u_old = Function(u, name="OldState")
m = [Function(W, name="Control_"+str(t)) for t in range(N)] # Create a control function for each timestep
v = TestFunction(V)

# Forward model
def run(u, u_old, m, annotate=True):
    u_old.vector()[:] = 0
    bc = DirichletBC(V, 0, "on_boundary")
    F = ((u - u_old)*v/timestep + inner(grad(u), grad(v)) + u**3*v - m[0]*v)*dx 
    # Perform N timesteps
    for t in range(N):
        F_t = replace(F, {m[0]: m[t]}) # Substitute the control function for the current timestep into the form
        solve(F_t == 0, u, bc, annotate=annotate)
        u_old.assign(u)

# Functional of interest
x = triangle.x
u_d = exp(-1/(1-x[0]*x[0])-1/(1-x[1]*x[1]))
u_dp = project(u_d, V, annotate = False)
J = Functional((0.5*inner(u-u_d, u-u_d))*dx*dt[FINISH_TIME])

# Optimise
run(u, u_old, m) # Create the tape
controls = [InitialConditionParameter(m[t]) for t in range(N)]
rf = ReducedFunctional(J, controls) 
m_opt = minimize(rf, options={"gtol": 1e-8}, bounds=([0]*N, [0.5]*N))
run(u, u_old, m_opt, False) # Compute the optimal state

# Plot
[plot(m_opt[t], title="Optimal control, t=%i"%t) for t in range(N)]
plot(u_dp, title="Desired state")
plot(u-u_d, title="Difference between optimal and desired state")