Optimizing State Value at Single Point in Time

dolfin-adjoint
1

I am interested in manipulating the source term of an ODE such that the state reaches a specified value at a specified time. For instance, in a minimal example, we would like to manipulate the in_weight function so that the state v is as close to 5 as possible at time point 6 (together with a norm regularization term):

J = (v(6)-5)**2*dx + 0.001*in_weight**2*dx

However, this gives me an error:

“UFLException: This integral is missing an integration domain.”

The error is thrown on that line.

I can optimize the source such that the ODE attains as close to that value as possible over the whole interval without issue:

J = (v-5)**2*dx + 0.001*in_weight**2*dx

The above line works fine.

I anticipate that the problem is that v(6) resolves to a floating point number, and the automatic adjoint chain is broken because there is no reference to the variable v. I do not know how to resolve this.

Here is a python script giving the error for a simple state equation:

from dolfin import *
from dolfin_adjoint import *
import numpy as np
import math
import time
import moola


def bd_func(shift):
        def boundary(x):
                return x[0] - shift < DOLFIN_EPS
        return boundary

V0 = Constant(0.0)

# Set up function space
mesh = IntervalMesh(1000, 0, 2*math.pi)
W = FunctionSpace(mesh,'CG',1)
bc = DirichletBC(W, V0, bd_func(0.0))
v = Function(W)
u = TestFunction(W)
in_weight = interpolate(Expression('0.0', degree=2), W)

weak_form  =  v.dx(0)*u*dx -v*u*dx  - in_weight*u*dx
    
Jac = derivative(weak_form, v, TrialFunction(W))
    
solve(weak_form==0,v, J=Jac,bcs=bc)

# Optimization at t = 6 (does not work)
J = (v(6)-5)**2*dx + 0.001*in_weight**2*dx
# Optimization over the whole interval (does work)
J = (v-5)**2*dx + 0.001*in_weight**2*dx

J=assemble(J)

rf = ReducedFunctional(J, Control(in_weight))

problem = MoolaOptimizationProblem(rf)
f_moola = moola.DolfinPrimalVector(in_weight)
solver = moola.BFGS(problem, f_moola)
sol = solver.solve()
2

Dear @nwycoff,

The reason that J = (v(6)-5)**2*dx + 0.001*in_weight**2*dx doesn’t work, is that v(6) and 5 are floats and therefore has no integration domain.

Recently, we added support of point evaluation to pyadjoint. For your example,
J = (v(numpy.array([6]))-5)**2 + assemble(0.001*in_weight**2*dx)
Does the trick.
To visualize the solution:

import matplotlib.pyplot as plt
plot(v.block_variable.saved_output)
plt.grid()
print("Actual value: ", v.block_variable.saved_output(6))
plt.savefig("v.png")