I’m trying out different changes in terms of optimization (functional, penalization, etc.), but don’t know what can cause the problem in terms of fenics. Do you have any thoughts about what’s wrong or in which direction should I go?
This does not have anything to do with the subdomain marking. Changing the problem still gives a nan-derivative.
from fenics import *
from fenics_adjoint import *
from collections import OrderedDict
import numpy as np
import itertools
set_log_active(False)
nu = Constant(1e-5)
data = Constant(2)
alpha = Constant(1e-3)
c_max = 30
c_min = 0
nx = ny = 20
mesh = UnitSquareMesh(nx, ny)
V = FunctionSpace(mesh, "CG", 1)
dt = Constant(0.1)
T = float(dt)
e = 2.71828
# Penalization parameters
u_max = 2.05
gamma = 100
tetha = 0.0016
# Chance contraint approximation parameters
m_1 = 1.0
m_2 = 0.4
tau = 1e-03
alpha = 0.9
class Obstacle(SubDomain):
def inside(self, x, on_boundary):
return (between(x[1], (0.3, 0.7)) and between(x[0], (0.3, 0.7)))
obstacle = Obstacle()
domains = MeshFunction("size_t", mesh, mesh.topology().dim())
domains.set_all(0)
#obstacle.mark(domains, 1)
#dx = dx(subdomain_data = domains)
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 + nu*inner(grad(u), grad(v)) - f*v)*(dx(0) + dx(1))
F = ((u - u_0)/dt*v + nu*inner(grad(u), grad(v)) - f*v)*dx
a, L = lhs(F), rhs(F)
bc = DirichletBC(V, 0, "on_boundary")
t = float(dt)
Theta = - 1 + alpha + (1+m_1*tau)/(1+m_2*tau*e**(1/tau*(-u_0 + u_max)))
j = 0
while t <= T:
f.assign(ctrls[t])
data.t = t
d.assign(interpolate(data, V))
solve(a == L, u_0, bc)
# Implement a trapezoidal rule
if t > T - float(dt):
weight = 0.5
else:
weight = 1
# j += (weight*float(dt)*assemble((u_0 - d)**2*dx)
# + gamma*0.125*float(dt)*assemble((Theta + ((Theta)**2 + tetha)**0.5)**2*dx(1)))
j += (weight*float(dt)*assemble((u_0 - d)**2*dx)
+ gamma*0.125*float(dt)*assemble((Theta + ((Theta)**2 + tetha)**0.5)**2*dx))
t += float(dt)
return u_0, d, j
u, d, j = solve_heat(ctrls)
# alpha = Constant(1e-3)
# regularisation = alpha/2*sum([1/dt*(fb-fa)**2*dx 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()]
rf = ReducedFunctional(J, m)
print(rf.derivative()[0].vector().get_local())
Please reduce your problem by removing parameters and see if the issue persists.
1 Like
You were absolutely right about subdomains, thank you!
I’ve tried to reduce the problem using Max function for penalization of Chance Constraints and replacing other parameters with their values:
from fenics import *
from fenics_adjoint import *
from collections import OrderedDict
import numpy as np
import itertools
from ufl.operators import Max
from ufl.operators import exp
set_log_active(False)
nu = Constant(1e-5)
data = Constant(2)
alpha = Constant(1e-3)
c_max = 30
c_min = 0
nx = ny = 20
mesh = UnitSquareMesh(nx, ny)
V = FunctionSpace(mesh, "CG", 1)
dt = Constant(0.1)
T = 2
# Penalization parameters
u_max = 2.05
gamma = 100
# Chance contraint approximation parameters
# m_1 = 1.0
# m_2 = 0.4
# tau = 1e-03
# alpha_CC = 0.9
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 + nu*inner(grad(u), grad(v)) - f*v)*dx
a, L = lhs(F), rhs(F)
bc = DirichletBC(V, 0, "on_boundary")
t = float(dt)
# Theta = - 1 + alpha_CC + (1+m_1*tau)/(1+m_2*tau*e**(1/tau*(-u_0 + u_max)))
# Theta = 0.1 + 1.001/(1+0.0004*exp(1000*(-u_0 + u_max)))
j = (0.5*float(dt)*assemble((u_0 - d)**2*dx)
+ 100*float(dt)*assemble((Max(0,(0.1 + 1.001/(1+0.0004*exp(1000*(-u_0 + u_max)))))**2*dx)))
while t <= T:
f.assign(ctrls[t])
data.t = t
d.assign(interpolate(data, V))
solve(a == L, u_0, bc)
# Implement a trapezoidal rule
if t > T - float(dt):
weight = 0.5
else:
weight = 1
j += (weight*float(dt)*assemble((u_0 - d)**2* dx)
+ 100*float(dt)*assemble((Max(0,(0.1 + 1.001/(1+0.0004*exp(1000*(-u_0 + u_max)))))**2*dx)))
t += float(dt)
return u_0, d, j
u, d, j = solve_heat(ctrls)
alpha = Constant(1e-3)
regularisation = alpha/2*sum([1/dt*(fb-fa)**2*dx 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()]
lower = [interpolate(Constant(c_min), V) for i in m]
upper = [interpolate(Constant(c_max), V) for i in m]
rf = ReducedFunctional(J, m)
print(rf.derivative()[0].vector().get_local())
But it doesn’t help.
I think exp(1/tau) must be the sorce of problem. Because this is the only thing, that caueses Overflow Exception by tau <= 1e03 computing it separetly:
exp(1000)
---------------------------------------------------------------------------
OverflowError Traceback (most recent call last)
<ipython-input-3-263c540ec676> in <module>
----> 1 exp(1000)
/usr/local/lib/python3.6/dist-packages/ufl/operators.py in exp(f)
608 def exp(f):
609 "UFL operator: Take the exponential of *f*."
--> 610 return _mathfunction(f, Exp)
611
612
/usr/local/lib/python3.6/dist-packages/ufl/operators.py in _mathfunction(f, cls)
593 def _mathfunction(f, cls):
594 f = as_ufl(f)
--> 595 r = cls(f)
596 if isinstance(r, (RealValue, Zero, int, float)):
597 return float(r)
/usr/local/lib/python3.6/dist-packages/ufl/mathfunctions.py in __new__(cls, argument)
108 def __new__(cls, argument):
109 if isinstance(argument, (RealValue, Zero)):
--> 110 return FloatValue(math.exp(float(argument)))
111 if isinstance(argument, (ComplexValue)):
112 return ComplexValue(cmath.exp(complex(argument)))
OverflowError: math range error
and it’s fully correlated with zero gradient (by the same tau values).
I think the following can handle this issue:
Function \Theta(\tau, \xi) = \frac{1+m_1*\tau}{1+m_2*\tau e^{\frac{1}{\tau}(u(\xi)-u_{max})}} \space (1) takes values near to zero if m_2 \tau e^{\frac{1}{\tau}(u(\xi)-u_{max})} \space (2) is high enough. Overflow Exception occurs also if (2) is too high, and I can basically set \Theta(\tau, \xi) to zero in this case.
exp(710) produces overflow exception. So if the condition m_2 \tau e^{\frac{1}{\tau}(u(\xi)-u_{max})} \geqslant e^{710} is sutisfied then \Theta(\tau, \xi) := 0 else compute Theta normally: \Theta(\tau, \xi) = (1)
But how can I implement this condition in terms of Dolfin adjoint, where Theta ( \Theta(\tau, \xi)) operates with u_0, and u_0 is a solution of PDE in 2D and ? I.e. Theta shouldn’t be just a constant value, but should be more like u_0.
You could use an UFL conditional:
from fenics_adjoint import *
from collections import OrderedDict
import itertools
nx = ny = 20
mesh = UnitSquareMesh(nx, ny)
V = FunctionSpace(mesh, "CG", 1)
dt = Constant(0.1)
T = float(dt)
ctrls = OrderedDict()
t = float(dt)
while t <= T:
ctrls[t] = Function(V)
t += float(dt)
def solve_heat(ctrls):
u_max = 2.05
u = TrialFunction(V)
v = TestFunction(V)
f = Function(V, name="source")
u_0 = Function(V, name="solution")
d = Function(V, name="data")
nu = Constant(1e-5)
F = ( (u - u_0)/dt*v + nu*inner(grad(u), grad(v)) - f*v)*dx
a, L = lhs(F), rhs(F)
bc = DirichletBC(V, 0, "on_boundary")
t = float(dt)
f.assign(ctrls[t])
d.assign(interpolate(Constant(2), V))
solve(a == L, u_0, bc)
# This works
j = (0.5*float(dt)*assemble((u_0 - d)**2*dx)
+ 100*float(dt)*assemble((0.1 + 1.001/(1+0.0004*exp(conditional(1000*(-u_0 + u_max)<500,1000*(-u_0 + u_max),500))))**2*dx))
# this does not work
# j = (0.5*float(dt)*assemble((u_0 - d)**2*dx)
# + 100*float(dt)*assemble((0.1 + 1.001/(1+0.0004*exp(1000*(-u_0 + u_max))))**2*dx))
return u_0, d, j
u, d, j = solve_heat(ctrls)
J = j
m = [Control(c) for c in ctrls.values()]
rf = ReducedFunctional(J, m)
print(rf.derivative()[0].vector().get_local())
Note from this example that a minimal example can be way smaller than what you created. You should think about this next time, as it would save people wanting to help you a considerable amount of time.
It works perfectly! Thank you so much for your support!
I will try to make my MWEs smaller.
I have one more question:
Do I understand correctly that this condition is imposed on every cell of the domain? Accordingly, condition output has different values in different cells depending on u_0 in these cells, but not constant value for the whole domain.
The conditional is evaluated at every quadrature point of the integral.
I suggest you play around with without solving optimization problems to get familiar with its behavior, i.e.
assemble(conditional(x[0]<0.5, 2*x[0], -x[1])*dx))
and consider what analytical output you expect from such an operation.
1 Like