R2015b PDE toolbox - Boundary conditions and solution of PDEs (PDE-constrained optimization Optimal Heating problem)

Hello, I am currently involved with PDE-constrained optimization in MATLAB R2015b.

I am trying to find a temperature control, u, on the boundary of a square [0:1] domain such that the temperature, y, in the domain comes as close to a desired temperature distribution, z_d, as possible without the control being unfeasible due to some control cost alpha.

The problem is to minimize the following cost function (I indicate the L2 norm in MATLAB sense, but I am aware that there is no norm(~,L2) in MATLAB). It is not important to my issue.

          J(y,u) = 1/2 * norm(y-z,L2)^2 + alpha/2 * norm(u,L2)^2
subject to:

          -Laplace y = 0                  in Omega
          dy/dn + rho*y = rho*u           on Gamma

By numerically estimating the gradient with gradJ = p+nu*u Where p is determined by solution of the following PDE system:

          -Laplace p = y-z_d              in Omega
           dp/dn = alpha*u                a.e. on Gamma
          -p*rho = alpha*u                on the rest of Gamma

For simplicity I want to solve the problem with the Dirichlet condition for p on the entirety of Gamma (-p*rho = alpha*u)

Note that rho, u, y, z_d are all functions of x = (x1, x2). Any rho and z_d can be chosen. I just want an example of how to do it and the choice of these functions aren't important.

I have tried different approaches and seem to like the approach of

model = createpde(1)
geometryFromEdges(model,geo)
msh = generateMesh(model)
[P, E, T] = meshToPet(msh)
applyBoundaryCondition(model,'Edge',...

However I seem to have made a mistake in my approach. So far my approach is:

• Define model for y and p, modely and modelp
• Generate the same mesh and so on for both models
• Generate initial guesses for u, p and y (zero vectors so far, u the length of the amount of boundary points, p and y the amount of points in the mesh)
• Apply boundary conditions to the respective models,
• Define the coefficients Cy, Ay and Fy for assempde(modely,Cy,Ay,Fy)
• Define the coefficients Cp, Ap and Fp for assempde(modelp,Cp,Ap,Fp)
• Solve for y
• Solve for p (here is where it goes wrong so far, I don't get any further)
• Invoke some kind of linesearch, for example sufficient decrease, and update u and boundary conditions
• Keep solving until tolerance level is satisfied

My questions are:

How does one define the neumann condition for y? dy/dn + y*rho = u*rho

I have tried

     applyBoundaryCondition(modely,'Edge',[1:modely.Geometry.NumEdges], 'g',ygfun ,'q',yqfun,'Vectorized','on');

Where ygfun = @(x1,x2) rho*intrpu and yqfun = @(x1,x2) rho(x1,x2)

Currently I have defined the neumann condition above 4 times, 1 on each edge and it seems to work however I cannot really check that yet until the p system at least runs.

I am having trouble defining the coefficients for the p system.

If you have any idea on how to define Cp, Ap and Fp that would be really helpful.

From the Poisson equation for p I get that Fp = y-z_d that is, a vector with the same length as the amount of mesh nodes.

But is this correct? Do I need to interpolate to the middle of each triangle? I don't really know at this point.

I have so much different, not-working code that I am afraid to include any of it here without a request.

[John BG]

This MATLAB CENTRAL contribution

http://uk.mathworks.com/examples/pde/community/6399-heating?s_tid=srchtitle

Seems to solve a heat transfer similar to the one you are on. Perhaps modifying the CAD model may help