Optimal control of heat source problem

dolfin-adjoint
1

Hi,
I try to minimize of the cost function with regulatization term

min_{u,m} 1/2 || u(T)-u_d ||^2 +\alpha/2 || m ||^2
subject to \frac{\partial u} {\partial t } - \nu*\Delta u = m
u=0 at time t=0 and u=0 on the boundary

I run the code and get the following error

    F = ((u-u_old)*v/dt + nu*inner(grad(u), grad(v)) - m[0]*v)*dx
Traceback (most recent call last):

  File "<ipython-input-66-4d4ef0e38cd7>", line 1, in <module>
    F = ((u-u_old)*v/dt + nu*inner(grad(u), grad(v)) - m[0]*v)*dx

  File "/usr/lib/python3/dist-packages/ufl/measure.py", line 437, in __rmul__
    domains = extract_domains(integrand)

  File "/usr/lib/python3/dist-packages/ufl/domain.py", line 351, in extract_domains
    return sorted(join_domains(domainlist))

TypeError: '<' not supported between instances of 'Mesh' and 'Mesh'

also get error for dt[FINISH_TIME]. It does not work neither. How can I make it work and find the plots.
the version of python is 3.6.9 and install on ubuntu.

The Code is

from fenics import *
import numpy as np
import matplotlib.pyplot as plt
from fenics_adjoint import *
from ufl import replace



N =30
dt = 0.2
T=1.0
nx = ny = 30



mesh = UnitSquareMesh(nx,ny)

V = FunctionSpace(mesh, "CG", 1)

W = FunctionSpace(mesh, "DG", 0)

u = Function (V, name="State")

u_old = Function(V, name="OldState")

x = SpatialCoordinate(mesh)

u_d =Expression("16*x[0]*(x[0]-1)*x[1]*(x[1]-1)*sin(pi*t)", t=0, degree=4)

nu= interpolate(Expression("1e-5",degree=2),V)

m = [Function(W, name="Control_"+str(t)) for t in range(N)]

v = TestFunction(V)

def run(m,annotate=False):
    u_old.assign(u_d)
    bc = DirichletBC(V, 0, "on_boundary")
    F = ((u-u_old)*v/dt + nu*inner(grad(u), grad(v)) - m[0]*v)*dx
    for t in range(N):
        F_t = replace(F,{m[0]:m[t]})
        solve(F_t == 0, u, bc, annotate=annotate)
        u_old.assign(u)



u_de = Expression("16*x[0]*(x[0]-1)*x[1]*(x[1]-1)*sin(pi*t)", t=1.0, degree=4)     

u_d = interpolate(u_de, V)

alpha = Constant(1e-7)

u_old=run(m,annotate=True)

J = assemble((0.5*inner(u_old-u_d,u_old -u_d))*dx*dt[FINISH_TIME] +0.5*alpha*sum(inner(m[i],m[i])*dx for i in range(N)))
 
controls = [Control(m[t]) for t in range(N)]

rf = ReducedFunctional(J, controls)

m_opt = minimize(rf, options={"disp": True})

u_opt = run(m_opt)

plot(u_d, title="optimised solution")

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

er = errornorm(u,u_d)

print(er)
2

Please revise your example, making sure it is runnable (Test it by copy-pasting your code into an empty file and run it). For instance nu is defined before the function-space, and the variable boundary is not defined.

3

So many thanks for your quick reply. I changed the code with the help of your advises. get the new error as following
Traceback (most recent call last):

  File "<ipython-input-68-0276cd086361>", line 1, in <module>
    runfile('/home/elif/Desktop/heatsource1.py', wdir='/home/elif/Desktop')

  File "/usr/lib/python3/dist-packages/spyder/utils/site/sitecustomize.py", line 705, in runfile
    execfile(filename, namespace)

  File "/usr/lib/python3/dist-packages/spyder/utils/site/sitecustomize.py", line 102, in execfile
    exec(compile(f.read(), filename, 'exec'), namespace)

  File "/home/elif/Desktop/heatsource1.py", line 63, in <module>
    J = assemble((0.5*inner(u_old-u_d,u_old -u_d))*dx*dt[FINISH_TIME] +0.5*alpha*sum(inner(m[i],m[i])*dx for i in range(N)))

TypeError: unsupported operand type(s) for -: 'NoneType' and 'Function'
4

That is because your run function doesn’t return anything. You have similar errors when trying to create sums of forms, and accessing dt[END_TIME], when dt is a float.
I’ve attached my modified and working version of the code. Please take more care when creating your minimal example next time, and spend some time trying to analyze your error. For instance use IPython.embed to enter interactive sessions inside your code.

from fenics import *
import numpy as np
import matplotlib.pyplot as plt
from fenics_adjoint import *
from ufl import replace

N =30
dt = 0.2
T=1.0
nx = ny = 30
mesh = UnitSquareMesh(nx,ny)
V = FunctionSpace(mesh, "CG", 1)
nu= interpolate(Expression("1e-5",degree=2),V)
W = FunctionSpace(mesh, "DG", 0)
u = Function (V, name="State")
u_old = Function(V, name="OldState")
x = SpatialCoordinate(mesh)
u_d =Expression("16*x[0]*(x[0]-1)*x[1]*(x[1]-1)*sin(pi*t)", t=0, degree=4)
bc = DirichletBC(V, 0, "on_boundary")
m = [Function(W, name="Control_"+str(t)) for t in range(N)]
v = TestFunction(V)

def run(m,annotate=False):
    u_old.assign(u_d)
    bc = DirichletBC(V, 0, "on_boundary")
    F = ((u-u_old)*v/dt + nu*inner(grad(u), grad(v)) - m[0]*v)*dx
    for t in range(N):
        F_t = replace(F,{m[0]:m[t]})
        solve(F_t == 0, u, bc, annotate=annotate)
        u_old.assign(u)
    return u_old


u_de = Expression("16*x[0]*(x[0]-1)*x[1]*(x[1]-1)*sin(pi*t)", t=1.0, degree=4)

u_d = interpolate(u_de, V)

alpha = Constant(1e-7)

u_old=run(m,annotate=True)

J = dt*0.5*inner(u_old-u_d,u_old -u_d)*dx
for i in range(N):
    J+= 0.5*alpha*inner(m[i],m[i])*dx
J = assemble(J)

controls = [Control(m[t]) for t in range(N)]

rf = ReducedFunctional(J, controls)

m_opt = minimize(rf, options={"disp": True})

u_opt = run(m_opt)

plot(u_d, title="optimised solution")

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

er = errornorm(u,u_d)

print(er)
5

Dokken, I really thank you!! You are so great!

6

I want to plot the m_opt for every time step but gets nothing, can you give an advice? btw is it necesarry to use dt in the following functional for my minimization problem?
If u run the code, u see that the error u_opt-u_d is so so small, is it possible?

J = dt*0.5*inner(u_old-u_d,u_old -u_d)*dx

from fenics import *
    import numpy as np
    import matplotlib.pyplot as plt
    from fenics_adjoint import *
    from ufl import replace
    from collections import OrderedDict

    T = 1
    dt = 0.1
    N = int(T/dt)
    m = OrderedDict()

    u_0 =Expression("16*x[0]*(x[0]-1)*x[1]*(x[1]-1)*sin(pi*t)", t=0, degree=4)

    mesh = UnitSquareMesh(10,10)

    V = FunctionSpace(mesh, "CG", 1)

    nu= interpolate(Expression("1e-5",degree=2),V)



    for t in range(N):
        m[t]=Function(V)
       

    def run(m,annotate=False):
        u = Function(V, name="State")
        v = TestFunction(V)
        u_old = Function(V, name="OldState")
        u_old.assign(u_0)
       
        F = ((u-u_old)*v/dt + inner(nu*grad(u), grad(v)) - m[0]*v)*dx 
        bc = DirichletBC(V, 0, "on_boundary")
        for t in range(N):
            F_t = replace(F,{m[0]:m[t]})
            solve(F_t == 0, u, bc, annotate=annotate)
            u_old.assign(u)
        return u_old

    u_0.t = T

    u_d = interpolate(u_0, V)

    alpha = Constant(1e-7)

    u_old=run(m, annotate=True)

    J = dt*0.5*inner(u_old-u_d,u_old -u_d)*dx
    for i in range(N):
        J += 0.5*alpha*inner(m[i],m[i])*dx
    J = assemble(J)

    controls = [Control(m[t]) for t in range(N)]

    rf = ReducedFunctional(J, controls)

    m_opt = minimize(rf, options={"disp": True})

    u_opt = run(m_opt)

    plot(u_opt, title="optimised solution")
    t=0
    for m in m_opt:
        plt.figure()
        t +=float(dt)
        plt.title("optimized source function at time= %g" %t)
        plot(m)
        plt.show()
        
      
    from matplotlib import pyplot, rc

    rc('text', usetex=True)

    x = [m((0.5, 0.5)) for m in opt_ctrls]

    pyplot.plot(x, label="$\\alpha={}$".format(float(alpha)))

    pyplot.ylim([-4, 4])

    pyplot.legend(loc=1)

     

    er = errornorm(u_opt,u_d)

    print("{:.20f}".format(er))
7

Also if I loosen the regularization coefficient alpha, I hope to see the control functions behave smoother but nothing happens. I am getting always the same m in every time period and for every alpha I tried. Why??

8

Please look through your example with care. Your desired profile u_d is equal to zero at t=1, and thus the optimal solution is 0, and all the controls will efficiently be zero.

9

Dokken, changing T didnot work, because I could only plot u_d but nothing else.

10

So what Did you change T to? You should look at u_D before posting any more questions about this. As your problem is not well formulated in mathematics, I cannot guide you any further. The optimal solution to the problem you posted above is the Zero vector, as u_D was equal to zero

11

sorry Dokken, I did not check u_D carefully, I try some couple of T and get similar u_D and u_opt to eachother. The error is not irrelivant as before. But again controls are zero.
I try the following example but try to change the cost function and regularization part.
http://www.dolfin-adjoint.org/en/latest/documentation/time-distributed-control/time-distributed-control.html
So many thanks for help Dokken.

12

Please produce a minimal working code that has a u_D\neq 0 that produces a minimal error and the control is 0 and producing an optimal solution.

13

Hello @dokken , thank you for your valuable contribution to @olric issue. I found it very useful for mu problem.

Please, how can I add the constraints on the control to the minimization routine?