Hello all,
I’m solving Time distributed control problem of heat equation, where I have some special term in functional. This term penalizes violation of the state contraint. For instance for state contraints u-u_{max} \geqslant 0 this term is \gamma(u-u_{max} + \sqrt{(u-u_{max})^2 + \theta})^2 . In deterministic case all works great. But I consider stochastic case where some parameters of PDE are random variables \xi, and in order to introduce state contraints I have to consider chance (CC) constraints of the form Pr\{u(\xi)-u_{max} \geqslant 0\} \geqslant \alpha. These CC I can approximate with the following function:
\Theta(\tau, \xi) = \frac{1+m_1*\tau}{1+m_2*\tau*e^{\frac{1}{\tau}(u(\xi)-u_{max})}} So that E[\Theta(\tau, \xi)] \leqslant 1 - \alpha . It all comes together in term I’m adding to functional: \gamma(\Theta(\tau) -(1 - \alpha) + \sqrt{(\Theta(\tau) - (1 - \alpha))^2 + \theta})^2 . (assuming I took only one sample from monte carlo simmulation). I won’t go into details any further, I just notice, that I have to solve this problem multiple times, each time decreasing \tau: \tau_{i+1} = \tau_i/10.
I have the following issue:
Expression e^{\frac{1}{\tau}(u-u_{max})} takes too high values for the python float data type with certain \tau. If I just try to calcualte this expression by \tau \geqslant 0.01 everything is alright, if \tau \leqslant 0.001 I have overflowerror, because if u-u_{max} =1 then it has e^{1000}, that can’t be computed. Talking about optimization, I can see the same: by \tau \geqslant 0.01 everything is alright, but by \tau \leqslant 0.001 at the first iteration projected gradient is zero and optimization doesn’t go further.
I tried to rescale the functional in order to increase the gradient, I tried of different L-BFGS_B tolerances (ftol, gtol), different step sizes (eps), I also tried TNC Algorithms, but projected gradient is always zero at the first iteration by \tau \leqslant 0.001.
I see two possible solutions:
- catching overflow exceptions, and setting \Theta to zero in this cases (because if overflow occurs, e^{\frac{1}{\tau}(u-u_{max})} \approx \infty \Rightarrow \Theta(\tau) = 0
- using decimal data type
But I don’t get any exceptions during optimization and I can’t get decimals work with fenics.
So I made minimal working example to ask about this issue. I just put control bounds and additional functional term into Time distributed control of heat equation demo, and suddenly it has worked by \tau \leqslant 0.001. (projected gradient is note zero at the first iteration). I adjusted optimization parameters a bit, and got to the right solution (near desired value and CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL).
But when I just run the code again, literally changing nothing, I don’t get the same result (algorithm converges too early CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH). I’ve reached good result multiple times, making small changes of the optimization parameters, but I could never reproduce it, optimization routine stops too early.
I have three questions:
- What could make my minimal working example so unstable? Usually I see convergence to the same solution, running code multiple times.
- What could be reason for projected gradient = 0 at the first iteration in my working code (I can laso provide it if needed) by \tau \leqslant 0.001.
- And how can one overcome these issues?
MWE:
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 = 2
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-02
alpha = 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 + (1+m_1*tau)/(1+m_2*tau*e**(1/tau*(-u_0 + u_max)))
j = (0.5*float(dt)*assemble((u_0 - d)**2*dx)
+ gamma*0.125*float(dt)*assemble((Theta + ((Theta)**2 + tetha)**0.5)**2*dx(1)))
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
Theta_1 = - 1 + alpha + (1+m_1*tau)/(1+m_2*tau*e**(1/tau*(-u_0 + u_max)))
j += (weight*float(dt)*assemble((u_0 - d)**2*dx)
+ gamma*0.125*float(dt)*assemble((Theta_1 + ((Theta_1)**2 + tetha)**0.5)**2*dx(1)))
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)
opt_ctrls = minimize(rf, method = "L-BFGS-B", options={"maxiter": 50,
'ftol': 2e-08,
'eps': 1e-02}, bounds = (lower,upper))
Good optimization results in MWE usually have functional value about 0.66 and look like this at each time step:
Thank you in advance for your time.