Error with Optimal control problem

General
1

Hi, I am trying to solve an optimal control problem using an example code I copied from here. However, I am getting errors that I don’t understand. Your help regarding this would be much appreciated. Thank you. Below is the code

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

Time parameters

N = 30
T = 2.0
num_steps = 5
dt = T/num_steps

Physical parameters

T_abs = 546.30 # Air absolute temperature in Kelvin
R = 8.31 # The ideal gas constant
w = 0.35 # Relative humidity of the air
P = 10e5; Air pressure

Summary

This text will be hidden

xa = (0.62198)/(0.62198 + w)
xw = w/(0.62198 + w)
rho_a = (28.9645P)(xa + 0.62198xw)/(RT_abs)
cp_a = 1 + 1.88*w

mesh = UnitSquareMesh(nx,ny)
V = FunctionSpace(mesh, “P”, 1)
W = FunctionSpace(mesh, “DG”, 0)
u = Function (V, name=“State”)
u_old = Function(V, name=“OldState”)
u_d = Function(V, name=“Desired state”)
x = SpatialCoordinate(mesh)
Q = [Function(W, name=“Control_”+str(t)) for t in range(N)]
v = TestFunction(V)
v_air = Function(W)

Define the initial solution and the air velocity vector valued

v_air = Constant((0.2, 0.2))
u_old = Constant(293.15)

Define boundary conditions

bc = DirichletBC(V, u_old, “on_boundary”)

As k_a depends upon the trial function T_a, it will be defined separately in the program as a function of T_a

def k_a (u):
return 0.02387 + 7.590e-5*u

def run(Q, annotate=False):
bc = DirichletBC(V, u_old, “on_boundary”)
F = uvdx + dt*(k_a(u)/(rho_acp_a))dot(v_air, grad(u))vdx + dt(k_a(u)/(rho_acp_a))dot(grad(u), grad(v))dx - dt(1/(rho_acp_a))Q[0]vdx - u_oldv*dx
for t in range(N):
F_t = replace(F,{Q[0]:Q[t]})
k_a(u_old)
solve(F_t == 0, u, bc, annotate=annotate)
u_old.assign(u)
return u_old

u_d = Constant(294.15)

alpha = Constant(1e-7)
gamma = Constant(1e-5)
y_min = Constant(20.0)
y_max = Constant(24.444)

u_old = run(Q,annotate=True)

J = dt0.5inner(u_old-u_d,u_old -u_d)dx + dt0.25gammainner(y_min - u_old + sqrt((y_min - u_old)**2), y_min - u_old + sqrt((y_min - u_old)**2))dx + dt0.25gammainner(u_old - y_max + sqrt((u_old - y_max)**2), u_old - y_max + sqrt((u_old - y_max)**2))dx
for i in range(N):
J+= 0.5alpha*inner(Q[i],Q[i])*dx
J = assemble(J)

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

rf = ReducedFunctional(J, controls)

Q_opt = minimize(rf,
bounds = [
[1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0],
[2000.0, 2000.0, 2000.0, 2000.0, 2000.0, 2000.0, 2000.0, 2000.0, 2000.0, 2000.0, 2000.0, 2000.0, 2000.0, 2000.0, 2000.0, 2000.0, 2000.0, 2000.0, 2000.0, 2000.0, 2000.0, 2000.0, 2000.0, 2000.0, 2000.0, 2000.0, 2000.0, 2000.0, 2000.0, 2000.0]],
options={“disp”: True})

u_opt = run(Q_opt)

fig1=plt.figure()

plt.subplot(1,2,1)

p = plot(u_d, title=“desired solution”)

cb = fig1.colorbar(p)

fig2=plt.figure()

plt.subplot(1,2,2)

p1 = plot(u_opt, title= “optmised solution”)

cb1 = fig2.colorbar(p)

#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)]
#fig=plt.figure()
#p = plot(u_d, title=“desired solution”)

fig3=plt.figure()

p3 = plot(m_opt[29])

cb3 = fig3.colorbar(p3)

#u_opt

er = errornorm(u_opt,u_d)

print(er)

2

Please make sure that you encapsulate your code with 3x`, i.e.

```python
from dolfin import *
# add other code
```

Please also state what error you are getting in a similar fashion:

```bash
Add error message here
```
3 
from dolfin import *
from fenics import *
import numpy as np
import matplotlib.pyplot as plt
from fenics_adjoint import *
from dolfin_adjoint import *
from ufl import replace

# Time parameters
N = 30
T = 2.0
num_steps = 100
dt = T/num_steps

# Physical parameters
T_abs = 546.30 # Air absolute temperature in Kelvin
R = 8.31 # The ideal gas constant
w = 0.35 # Relative humidity of the air
P = 10e5;
xa = (0.62198)/(0.62198 + w)
xw = w/(0.62198 + w)
rho_a = (28.9645*P)*(xa + 0.62198*xw)/(R*T_abs) # Air density
cp_a = 1 + 1.88*w # Heat capacity of air 

# The domains
nx = 30
ny = 20
mesh = UnitSquareMesh(nx,ny)
V = FunctionSpace(mesh, "P", 1)
W = FunctionSpace(mesh, "DG", 0)

# Functions definition
u = Function (V, name="State")
u_old = Function(V, name="OldState")
u_d = Function(V, name="Desired state")
x = SpatialCoordinate(mesh)
Q = [Function(W, name="Control_"+str(t)) for t in range(N)]
v = TestFunction(V)
v_air = Function(W)

# Define the initial solution and the air velocity vector valued
v_air = Constant((0.2, 0.2))
u_d = Expression("294.15", degree=1)
u_old = interpolate(u_d, V)
u_de = interpolate(u_d, V)  // The desired state

# Define boundary conditions
bc = DirichletBC(V, u_old, "on_boundary")

# As k_a depends upon the trial function T_a, it will be defined seperately in the program as a function of T_a
def k_a (u):
    return 0.02387 + 7.590e-5*u

def run(Q, annotate=False):
    bc = DirichletBC(V, u_old, "on_boundary")
    F = u*v*dx + dt*(k_a(u)/(rho_a*cp_a))*dot(v_air, grad(u))*v*dx + dt*(k_a(u)/(rho_a*cp_a))*dot(grad(u), grad(v))*dx - dt*(1/(rho_a*cp_a))*Q[0]*v*dx - u_old*v*dx
    for t in range(N):
        F_t = replace(F,{Q[0]:Q[t]})
        k_a(u_old)
        solve(F_t == 0, u, bc, annotate=annotate)
        u_old.assign(u)
    return u_old

alpha = Constant(1e-7)
gamma = Constant(1e-5)
u_min = Expression("293.15", degree=1)   // Lower bound of the state
u_max = Expression("297.594", degree=1)  // Upper bound of the state

u_old = run(Q,annotate=True)

J = dt*0.5*inner(u_old-u_de,u_old -u_de)*dx + dt*0.25*gamma*inner(y_min - u_old + sqrt((y_min - u_old)**2), y_min - u_old + sqrt((y_min - u_old)**2))*dx + dt*0.25*gamma*inner(u_old - y_max + sqrt((u_old - y_max)**2), u_old - y_max + sqrt((u_old - y_max)**2))*dx
for i in range(N):
    J+= 0.5*alpha*inner(Q[i],Q[i])*dx
J = assemble(J)

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

rf = ReducedFunctional(J, controls)

Q_opt = minimize(rf, 
                 bounds = [
                [1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0],
                [2000.0, 2000.0, 2000.0, 2000.0, 2000.0, 2000.0, 2000.0, 2000.0, 2000.0, 2000.0, 2000.0, 2000.0, 2000.0, 2000.0, 2000.0, 2000.0, 2000.0, 2000.0, 2000.0, 2000.0, 2000.0, 2000.0, 2000.0, 2000.0, 2000.0, 2000.0, 2000.0, 2000.0, 2000.0, 2000.0]], 
                options={"disp": True})

u_opt = run(Q_opt)

fig=plt.figure()
p = plot(u_opt, title="Optimised solution")
cb = fig.colorbar(p)

#plot(u_de, title="Desired solution")

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

er = errornorm(u_opt ,u_de)
print(er)

plt.show()

Hi dokken, many thanks for your reply. I was able to fix the previous error. Now I have another concern about the result I am getting. While plotting the optimised state, the values are not varying that much and it’s almost constant. Moreover the error is satisfactory (err = 0.06405098730435595). I don’t what can be the problem in this case. Could you please help me to understand it a bit clearer? Thank you

4

Sorry! I meant the error is not satisfactory. Thanks