How to set state constraints in PDE constrained optimal control

dolfin-adjoint
1

Hello everyone,

I’m trying to set temperature constraints in Optimal control of heat equation. As I know, one way to do it is penalization technique. So I tried Moreau-Yosida regularization and didn’t get the goal. (state is steal not constrained).

I was wondering, what else can I do to set these constraints or what should be the right way?

Here is my initial code without attemps to use Moreau-Yosida regularization.

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

data = Expression("16*x[0]*(x[0]-1)*x[1]*(x[1]-1)*sin(pi*t)", t=0, degree=4)
nu = Constant(1e-5)

nx = ny = 20
mesh = UnitSquareMesh(nx, ny)
V = FunctionSpace(mesh, "CG", 1)

class Left(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0.0)
class Right(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 1.0)
class Bottom(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 0.0)
class Top(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 1.0)
    
left = Left()
top = Top()
right = Right()
bottom = Bottom()
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundaries.set_all(0)
left.mark(boundaries, 1)
top.mark(boundaries, 2)
right.mark(boundaries, 3)
bottom.mark(boundaries, 4)
g1 = Constant("0.01")
g2 = Constant("0.01")
g3 = Constant("0.01")
g4 = Constant("0.01")
ds = ds(subdomain_data = boundaries)

class Obstacle(SubDomain):
    def inside(self, x, on_boundary):
        return (between(x[1], (0.4, 0.6)) and between(x[0], (0.4, 0.6)))
    
obstacle = Obstacle()
domains = MeshFunction("size_t", mesh, mesh.topology().dim())
domains.set_all(0)
obstacle.mark(domains, 1)
dx = dx(subdomain_data = domains)

dt = Constant(0.1)
T = 2
steps = 20

ctrls = OrderedDict()
t = float(dt)
while t <= T:
    ctrls[t] = Function(V)
    t += float(dt)

def solve_heat(ctrls):
    u = TrialFunction(V)
    v = TestFunction(V)
    f = Function(V, name="source")
    u_0 = Function(V, name="solution")
    d = Function(V, name="data")
    F = ((u - u_0)/dt*v*(dx(0) + dx(1)) + nu*inner(grad(u), grad(v))*dx(0) + nu*inner(grad(u), grad(v))*dx(1) - f*v*dx(1) - g1*v*ds(1) - g2*v*ds(2) - g3*v*ds(3) - g4*v*ds(4))
    a, L = lhs(F), rhs(F)
    t = float(dt)
    j = 0.5*float(dt)*assemble((u_0 - d)**2*(dx(0) + dx(1)))
    while t <= T:
        f.assign(ctrls[t])
        data.t = t
        d.assign(interpolate(data, V))
        solve(a == L, u_0)

        if t > T - float(dt):
           weight = 0.5
        else:
           weight = 1

        j += weight*float(dt)*assemble((u_0 - d)**2*(dx(0) + dx(1)))
        t += float(dt)

    return u_0, d, j

u, d, j = solve_heat(ctrls)

alpha = Constant(1e-2)
regularisation = alpha/2*sum([1/dt*(fb-fa)**2*(dx(0) + dx(1)) for fb, fa in 
                              zip(list(ctrls.values())[1:], list(ctrls.values())[:-1])])

J = j + assemble(regularisation)
m = [Control(c) for c in ctrls.values()]

c_max = 1.5
c_min = -1.5
lower = [interpolate(Constant(c_min), V) for i in range(steps-1)]
upper = [interpolate(Constant(c_max), V) for i in range(steps-1)]

rf = ReducedFunctional(J, m)
opt_ctrls = minimize(rf, options={"maxiter": 50}, bounds = (lower,upper))

Appreciate any help.

2

Regularization is usually the way to go as you are working with a temporal finite difference scheme.