Just to make sure I'm getting the correct plots for the optimals

Asked by ahmad ali

Hallo All
I want to have some graphs for the optimal control and the optimal solution for my problem. Thanks to the code
in https://gist.github.com/funsim/5769594 I got some plots.

here is my code

from dolfin import *
from dolfin_adjoint import *

# Setup
n = 40
mesh = UnitCircleMesh(n)
V = FunctionSpace(mesh, "CG", 1)
###W = FunctionSpace(mesh, "DG", degree = 0)

u = Function(V,name="State")
m = Function(V, name="Control")
v = TestFunction(V)

# Run the forward model once to create the simulation record
F = ( inner(grad(u), grad(v)) + u*v - m*v )*dx
###bc = DirichletBC(V, 0.0, "on_boundary")
solve(F == 0, u)

def max(x, y):
  import ufl
  x = triangle.x[0]
  return ufl.conditional(gt(x, y), x, y)

# the alpha and gamma coefficients which appear in the cost functional
alpha=Constant(1.0)
gamma=Constant(100000.0)

# The given u_0 , m_0 and b
x = triangle.x
u_0 = 4+1/pi -(1/(4*pi))*(x[0]*x[0]+x[1]*x[1])+(1/(2*pi))*ln(sqrt(x[0]*x[0]+x[1]*x[1]))*(1/ln(10)) #(exp(-1/(1-x[0]*x[0])-1/(1-x[1]*x[1]))
m_0 = 4 +(1/(4*1*pi))*(x[0]*x[0]+x[1]*x[1])-(1/(2*1*pi))*ln(sqrt(x[0]*x[0]+x[1]*x[1]))*(1/ln(10))
b= (x[0]*x[0]+x[1]*x[1])+4

# The cost functional

J = Functional(( 0.5*inner(u-u_0, u-u_0) + alpha*0.5*inner(m-m_0,m-m_0) + 0.5*(1/gamma)*inner(max(0,gamma*(u-b)),max(0,gamma*(u-b))) )*dx)

# Run the optimisation
reduced_functional = ReducedFunctional(J, InitialConditionParameter(m, value=m))
# Make sure you have scipy >= 0.11 installed
m_opt = minimize(reduced_functional, method = "L-BFGS-B",
                 tol=2e-08, options = {"disp": True})

# desire solution u_0
proj_udes = project(u_0,V)
plot(u_0, mesh, interactive=True, title="u_0")
# Actual solution
solve(F == 0, u)
plot(u, interactive=True, title="Optimal State")
# L1 error
proj_error = project(u-u_0,V)
plot(proj_error, interactive=True, title="L1-Error")
# Optimal forcing
plot(m_opt, interactive=True, title="Optimal Control")

regardless of the problem or its parameters, I want to make sure if this code is really giving plots for the optimal control and optimal state.
  I'm assuming if one wants to plot the optimal state , then one needs to slove a new bilinear form which has the optimal control.

Just I was testing the software on some problems with known solutions and I found a slight diferece.
so, I wanted to make sure if my code is really written in a way that gives the plots of the optimals.

Question information

Language:
English Edit question
Status:
Solved
For:
dolfin-adjoint Edit question
Assignee:
No assignee Edit question
Solved by:
Simon Funke
Solved:
Last query:
Last reply:
Revision history for this message
Best Simon Funke (simon-funke) said :
#1

Hi Ahmed,

Looking good!

A few suggestions:

* Try "Newton-CG" as your optimisation algorithm. You will see it converges quicker.
* The optimal control has visible noise in its solution. This is caused by the fact that scipy.optimise does not know anything about the function spaces that your optimisation parameterse live in - and uses the Euclidian inner product (instead of the proper L2 inner product) during its computations. You can force dolfin-adjoint to make a conversion step for scipy by adding the parameter "in_euclidian_space=True" to "minimize". However, you will see that the current implementation is really slow. Since the Newton-CG algorithm does not suffer from this, so I suggest that you just stick with "Newton-CG" for now until I have implemented alternatives for the scipy optimisation algorithms.

> I'm assuming if one wants to plot the optimal state , then one needs to slove a new bilinear form which has the optimal control.
Yes.

Simon

Revision history for this message
ahmad ali (cfcf1111) said :
#2

Thanks Simon Funke, that solved my question.