Lesson 5

Lesson 5 : An Optimization Problem

The present lesson demonstrates how to use an optimization problem solver. We solve the same problem as in Lesson 2 as an optimization problem where Dirichlet boundary conditions are considered as equality constraints. The optimization algorithm is the Truncated Newton algorithm described in function OptimTN. This is not the best method to solve a Laplace equation, but our purpose here is to learn how to use this class.

The Finite Element Code

The main program

Here is a description of the source code. The source file can be found in Example 5 in the examples directory.

We start, as usual, by including required headers and naming the appropriate namespace. We furthermore include the header file of the optimization definition class called Opt.
```
#include "OFELI.h"
#include "Opt.h"
#include "User.h"
using namespace OFELI;
```
The program will have as argument the name of the parameter data file :
```
int main(int argc, char *argv[])
{
```
We next declare a pointer to class Element that will be used later.
```
   Element *el;
```

We expand program arguments and declare an instance of class IPF for parameter file :

   if (argc <= 1) {
      cout << " Usage : ex5 <parameter_file>" << endl;
      exit(1)
   }
   IPF data(argv[1]);

The mesh data file is obtained from a member function of instance data.
```
   Mesh ms(data.MeshFile());
```
We introduce boundary conditions through a user defined class (See Lesson 4).
```
   User ud(ms);
```
Implementation of class User will be given later.
n is the number of degrees of freedom.
```
   int n = ms.getNbDOF();
```
We declare some vectors. x is the solution vector, vectors low and up will contain for each degree of freedom upper and lower bound respectively to prescribe and pivot is a vector that will contains (after optimization) flags to indicate which constraint is reached.
```
   Vect<double> x(n),low(n),up(n);
   Vect<int> pivot(n);
```
Dirichlet boundary conditions are taken into account as before :
```
   Vect<double> bc(n);
   ud.setDBC(bc);
```
We start now to properly define the optimization problem. For this, we declare an instance of a class Opt that is defined separately. The constructor of this class invokes the mesh instance and the instance ud that will be useful to transmit to the class body or boundary sources.
```
   Opt theOpt(ms,ud);
```
We also initialize the solution to zero :
```
   x = 0;
```
As said before, Dirichlet boundary conditions are considered here as equality constraints and are then incorporated into vectors low and up via a utility function called BCAsConstraint contained in the mentioned file OptimAux.h
```
   BCAsConstraint(ms,bc,up,low);
```
The function OptimTN can now be invoked to run the optimization algorithm.
```
   OptimTN<Opt>(theOpt,x,low,up,pivot,100,1.e-12,1);
```
The reader can refer to the function OptimTN to understand the meaning of each argument of the function.
If the algorithm has succeeded (which is the case in this example) we can output the solution and close the program.
```
   cout << "\nSolution :\n" << x;
}
```

A User defined class

We have now to implement class Opt that defines the problem to solve and provides to the objective function and its gradient. The class is defined in file Opt.h that we study here below.

We start by including files OFELI.h and Therm.h since we are going to solve a heat transfer problem. We also invoke the namespace OFELI :
```
#include "OFELI.h"
#include "Therm.h"
#include "User.h"
using namespace OFELI;
```
This class will have a constructor that acquires the mesh and the user data instance and stores pointers to these objects:
```
class Opt {

public:
   Opt(Mesh &ms, User &ud) { _ms = &ms; _ud = &ud; }
```
The other public member is the objective function.
```
   void Objective(Vect<double> &x, double &f, Vect<double> &g) {
        Element *el;
        f = 0.;
        g = 0.;
        for (_ms->topElement(); (el=_ms->getElement());) {
           DC2DT3 eq(el);
           Vect<double> ge(3);
           Vect<double> xe(el,x);
           f += eq.Energy(xe,*_ud);
           eq.EnerGrad(xe,*_ud,ge);
           g.Assembly(el,ge);
        }
   }
```
Note that the arguments of the objective function are necessarily the optimization variable vector, the objective function to compute and the gradient vector to compute. Here the class DC2DT3 provides necessary material for optimization purposes. Namely, member function Energy return the energy value and function EnerGrad calculates the element energy gradient that needs be assembled into the global vector g.
It remains to declare the mesh and user data pointers as private attributes of the class and end the class definition.
```
private:
   Mesh *_ms;
   User *_ud;
};
```

A test

We use here exactly the same mesh file as in the previous lesson. Of course, we obtain the same solution as the previous lesson. The output of the program is the following:



   NIT  NF   CG           F                  GTG

    0    1    0     6.00000000e+000     4.00000000e+001
    1    2    3     1.67852990e+000     4.09059524e+001
    2    3    5     1.50042385e+000     4.09106651e+001
    3    4    8     1.50001179e+000     4.09109024e+001
    4    5   10     1.50000062e+000     4.09109074e+001
    5    6   13     1.50000000e+000     4.09109074e+001
    6    7   15     1.50000000e+000     4.09109074e+001
    7    8   17     1.50000000e+000     4.09109074e+001
    8    9   19     1.50000000e+000     4.09109074e+001


Optimal Function Value = 1.5

Solution :
     1     1.00000000e+000
     2     7.49999986e-001
     3     4.99999972e-001
     4     2.49999994e-001
     5     0.00000000e+000
     6     1.00000000e+000
     7     7.49999988e-001
     8     4.99999982e-001
     9     2.49999978e-001
    10     0.00000000e+000
    11     1.00000000e+000
    12     7.49999989e-001
    13     4.99999982e-001
    14     2.49999978e-001
    15     0.00000000e+000
    16     1.00000000e+000
    17     7.49999986e-001
    18     4.99999972e-001
    19     2.49999994e-001
    20     0.00000000e+000