Optimization of Finite Element Models
with MathOptimizer and ModelMaker

János D. Pintér
Pintér Consulting Services, Inc., Halifax, NS, Canada 
jdpinter@hfx.eastlink.ca
and
Christopher J. Purcell
Defence Research and Development Canada-Atlantic, Dartmouth, NS, Canada
chris.purcell@drdc-rddc.gc.ca

Introduction

Quantitative decisions related to engineering, economic and scientific investigations are frequently made using optimization concepts and tools. The decision-maker or modeler typically wants to find the "absolutely best" decision which corresponds to the minimum or maximum of a suitable objective function and satisfies certain feasibility constraints. The objective function expresses overall system performance and the constraints originate from physical, technical, or economic considerations.

In this presentation we shall consider engineering applications of the continuous global optimization (GO) model stated as 

	min f(x)	   f: (Subscript[D, 0])->(R^1)	
	g(x)=0		g: (Subscript[D, 0])->(Subscript[(R^m), e])
	h(x)≤0		h: (Subscript[D, 0])->(Subscript[(R^m), i])
	xl≤x≤xu	  x, xl, xuϵ(R^n)
	f, g, h are continuous functions. 

All equality and inequality constraints are interpreted component-wise.

In its full generality, the GO model encompasses a large class of very difficult instances. Figure 1 below illustrates the potential difficulty of continuous optimization models. Although the feasible set is just a two-dimensional interval, the objective function is highly multimodal. Hence, purely local scope search methods—as a general rule,  even in much simpler cases—will not find the best solution.

Figure 1. A box-constrained global optimization model in two dimensions.

MathOptimizer

Mathematica has several built-in functions for solving certain types of optimization models: these include e.g., LinearProgramming (for linear programming), FindMinimum for unconstrained (nonlinear) convex optimization, and NMinimize for—in principle, arbitrary—constrained optimization.

MathOptimizer is a third-party Mathematica package that enables the global and local solution of the general class of optimization problems encompassed by the GO model form. The key features of this package are summarized below: for theoretical background, see [1, 2].

MathOptimizer currently consists of two core solver packages and a solver integrator package. One of the solver components serves to the approximate global optimization of an aggregated merit function on the given interval range [xl,xu]: the merit function includes the objective function f and penaltes for the (possibly) violated constraints g and h. This package is based on an efficient adaptive stochastic search method, combined with a statistical bounding procedure: the latter provides an estimate of the global optimum value, based on the global sampling results.

The second solver package is aimed at finding a vector x that satisfies the Karush-Kuhn-Tucker (local) optimality conditions. This solver option can be used for "precise" local optimization, based on a given initial solution that is either produced by a preceding global search phase, or is directly provided by the user.

The solver integrator package supports the individual or combined use of these two solver packages; the solvers can also be used separately.

The MathOptimizer User Guide [3] is a Mathematica notebook of over 70 (printed) pages that can be directly invoked through the on-line Help system. The User Guide includes installation and technical notes, and provides concise mathematical background information and modeling tips. In addition to basic usage description, it also discusses a number of test problems and applications, including the challenging GO model example shown in Figure 1. This example is based on Problem 4 of A Hundred-dollar, Hundred-digit Challenge, posted by Trefethen [4]: find the global minimum of the function

Exp[Sin[50*x]] + Sin[60*Exp[y]] + Sin[70*Sin[x]] + Sin[Sin[80*y]] - 
        Sin[10*(x + y)] + (x^2 + y^2)/4;

The solution found by MathOptimizer has an estimated optimum value of -3.306868647... this result is precise to at least ten digits. (The corresponding runtime is a fraction of a second on a P4 1.6 GHz processor-based machine that runs Mathematica 5.0 under Windows.)

ModelMaker

ModelMaker, a package for finite element modelling, has been described in previous Developer conferences [5, 6]. It builds parameterized models for the design of sensors, vehicles or structures. Parameters can specify geometry, material properties, loads, and meshing rules. The models are passed to external analysis engines for processing and results are imported for interpretation. The entire process can be coded as a single Mathematica function, which is convenient for design iterations. These could be done by trial and error, but also can readily be put under the control of an optimizer. MathOptimizer was developed, in part, to provide a state-of-art nonlinear (global) optimizer to use with ModelMaker.

The design optimization example presented below is from acoustical engineering; note that ModelMaker and MathOptimizer have much wider applicability. The optimization tools used here are "overkill" for this one-parameter problem, but we wished to provide a very simple illustration of how ModelMaker is interfaced with an optimizer.

Example: Design of a Tuning Fork

This example designs a tuning fork with a prescribed fundamental frequency, by varying its height. The fork has a square cross section. The figure below shows the finite element mesh - using 20 noded hexahedral elements.  We model the fork as an eigenvalue problem in vacuum, neglecting the small frequency shift due to acoustic radiation damping that a tuning fork vibrating in air will experience. The finite element analysis is done with a code called Mavart3D [7]. Data is passed to Mavart3D by writing to an external file, the code is then invoked with ModelMaker's Analyze command, and the results of analysis are imported back into the notebook.

If the ModelMaker package is in $AddOnsDirectory, you install it and the definitions that support the Mavart3D finite element code by executing the next cell.

\(\(SetDirectory[
        ToFileName[{$AddOnsDirectory, "Applications", "ModelMaker"}]];\)\)
Needs["ModelMaker`mm`"]; 
Get["ModelMaker`mavart3d`"];

Then you can run the following cell which builds the geometry of the fork using some initial estimates for the dimensions, meshes it with parabolic hex elements, writes out the model, runs the analysis to calculate the eigenfrequencies (in Hertz) of the tuning fork, and imports the results. The first non-rigid-body mode is 7, since the fork is unconstrained, so the first 6 trivial solutions are dropped.

\(\(\(Point[Delete]; Node[Delete]; Fixity[Delete]; \ 
    Material[Delete];  (*\ Dimensions\ in\ meters\ *) 

    Clear[h, r, w];  (*\ 
      the\ value\ of\ height\ h\ will\ be\ set\ later\ *) 

    r\  =  .015; \ \  (*\ radius\ of\ the\ fork\ *) 

    w\  = 0.009; \ \  (*\ width\ of\ tines, 
      fork\ has\ a\ square\ cross - section\ *) 

    Point[New, {{0, 0}, {0, h}, {0, h + w}}]; \  (*\ 
      define\ cross\ section\ in\ the\ x - y\ plane\ *) 

    Point[Translate, {1, 2, 3}, {w/2, 0}, 1];

    \(Point[Rotate, {5, 6}, {0, 0, \(-45\)\ Degree}, 
        Point[\(Point[Translate, Point[Last], {0, r}, 1]\)[\([1]\)]\ ], 
        2];\)

    \(Point[Translate, {10, 11}, {0, h}, 1];\)

    Curve[New, {1, 4}]; Curve[New, {2, 5}]; Curve[New, {3, 6}];

    Curve[New, {5, 8, 10}]; Curve[New, {6, 9, 11}]; Curve[New, {11, 10}]; 
    Curve[New, {13, 12}];

    Surface[New, {1, 2}]; Surface[New, {2, 3}]; Surface[New, {4, 5}]; 
    Surface[New, {6, 7}];

    \(Surface[Translate, SurfaceList, {0, 0, w}, 1];\)

    Solid[New, {1, 5}]; Solid[New, {2, 6}]; Solid[New, {3, 7}]; 
    Solid[New, {4, 8}];

    Node[Delete]; Mesh[2, SolidList, {3, 3, 1}];\ \ 

    Element[Mirror, "X"]; Node[Equivalence];

    Material[New, 3, {7700. , 0, 0, 0, 190.*^9, 0.3}, "Steel"]; \  (*\ 
      Density, \ Young' s\ modulus, \ 
      Poisson\ ratio\ *) 

h\  =  .1; \  (*\ 
      set\ the\ height\ before\ writing\ data\ file\ *) 

\ 
    WriteModel["mav.dat", 
      TITLE -> 
        StringJoin["tuning fork ", 
          ToString[Length[NodeList]], " Nodes "], 
      MACROS -> {"PRNT EIGE"}, RUNTASK -> {"EIGE"}, 
      NF -> 12];


    \(\(Analyze[]\) \); \)\(\ \)\( (*\ 
      Run\ analysis\ and\ retrieve\ results\ *) \)\)\)

TableForm[
      Transpose[Transpose[
            Drop[
ReadResult["mav.out", 
                ResultType -> "EigenValues"], 6]][[{1, 
              2}]]]]

Our inspired guess (inspired by looking at a real A440 tuning fork) has yielded a lowest natural resonance frequency of 422.5 Hz. Next we import the eigenvector for Mode 7, and plot an animation. It looks like what you would expect for the lowest mode of the tuning fork. Double click the graphic below to play the animation.

ReadResult["mav.out",ModeNumber->7];
Element[Plot,,DeformedElementForm->"Fill",ViewVertical->{0,1,0},
ViewPoint-> {1,1,1},EdgeWidth->0.001,
PlotRange-> 1.1 * Node[Maximum]{{-.5,.5},{0,1},{-.1,.1}},
DeformSize->5,Boxed->True,Axes->{True,True,False},AxesLabel->{"x","y","z"},
Deform->True, Undeform->False,ElementNodes->False,Frames->10];

Using Optimization to Tune the Fork

If we clear the value of the height h and list the coordinates of the first 10 nodes we see they contain both numbers and expressions involving h. This illustrates a key feature of ModelMaker - the finite element model can contain symbolic expressions which enable rapid design changes.

Clear[h]; Node[#][[3]]& /@ Range[10]

To make design iterations, we can assign values to the symbol h, export the model, run the analysis and examine the eigenfrequencies. The geometry and mesh do not have to be regenerated, which yields a considerable savings of time.

Next, with the symbolic model still in memory, we define the objective function that the optimizer will minimize, to find the length of the tuning fork for 440 Hz. To display the convergence, the objective function, as a side effect, creates a notebook with a plot of resonance frequency versus fork length, and updates the plot on each objective function evaluation. The objective function evaluation results are collected in the symbol evals. The most recent versions of FindMinimum and NMinimize, in Mathematica Version 5.0, have the options StepMonitor and EvaluationMonitor, and can use the routines Reap and Sow, to provide alternative tools for monitoring progress.

evals = {}; (* list of all objective function evaluations *)
target=440.0; (* the goal for the optimization - frequency in Hz of mode 7 *)

Clear[obj,h];  
obj[x_?NumberQ] := Module[{fr,nb}, (
      Notebook[{Cell["Output from Optimization", "Text"]}] >> "output.nb"; 
      nb = NotebookOpen[ToFileName[Directory[],"output.nb"], WindowSize -> {500, 400}]; 
      h=x;   (* insert the height of the tines into the model *)   
      WriteModel["mav.dat", TITLE -> StringJoin["tuning fork ",ToString[Length[NodeList]],
       " Nodes "], 
      MACROS -> {"PRNT EIGE"}, RUNTASK -> {"EIGE"}, NF -> 7]; 
      Analyze[]; 
      fr =Drop[ReadResult["mav.out",ResultType->"EigenValues"],6][[1,2]];
      AppendTo[evals, {h, fr}]; 
      ListPlot[evals, Frame -> True, RotateLabel -> False, FrameLabel -> {"height","Hz"}, 
      Axes -> False, 
      PlotStyle -> {PointSize[0.03], Hue[1]}, DefaultFont -> {"HelveticaBold",10}, 
      PlotRange->{{0,.2},{0,3 target}},GridLines -> {None,{target}},
      PlotLabel -> StringJoin[ToString[Length[evals]],
      "  height = ",ToString[PaddedForm[h, 4]],
      "     f1 = ", 
      ToString[PaddedForm[fr, 6]], " Hz"], 
      DisplayFunction -> (SelectionMove[nb,All,Notebook]; 
      NotebookWrite[nb, Cell[GraphicsData["PostScript", 
      DisplayString[#1]], "Graphics"], All] & )]; 
      Abs[target - fr]
    )];

Next we show that with this objective function, Mathematica's built-in function FindMinimum can find the height h that meets our design goal. The graphic below is displayed in a separate notebook, and shows the progress of the search for a design. The horizontal line is at an ordinate of 440 Hz. The result is that 30 objective function evaluations are sufficent to find a design that gets to within 0.1 Hz of our goal. The tines of the fork must be cut to a length of  0.097 meters. The Timing results shown are for a G4 processor at 867 MHz clock speed.

evals={}; (* initialize the list of results *)
Timing[
  FindMinimum[obj[Height],{Height,0.05,.15}] ]
Length[evals];

Next we reset the evaluation history, and use NMinimize to solve the problem. It needs 1723 objective functions evaluations to satisfy the design goal, using the default settings for all its options.

evals={};
Timing[NMinimize[{obj[Height],0.05<Height<.15},Height]]
Length[evals]

Finally we reset the evaluation history, load MathOptimizer, and repeat the above calculation using values for options that activate the global search mode, and deactivate the gradient based local search (which can't be used here). It returns a solution that is within 0.1 Hz of our goal while requiring 702 function evaluations. This objective function has only one local minimum in the search region, so the effort MathOptimizer expends in checking for local minima is unnecessary for this simple problem.

Needs["MathOptimizer`Optimize`"];
evals = {};
Timing[Optimize[obj[Height],{},{},{Height},{.05},{.08},{.15}, 
      GlobalSolverMode->1, LocalSolverMode->0]]
Length[evals]

Optimization of Expensive "Black Box" Models

MathOptimizer's current local solver is intended for objective functions with first derivatives that can be evaluated by Mathematica in symbolic form: this approach can't be used when the objective function is a black box like in the tuning fork example. MathOptimizer is also intended for objective functions which are low cost to evaluate, because the random search method it employs in global search mode typically requires a large number of function evaluations. The next upgrade to MathOptimizer, scheduled for fall 2004, will add a capability to deal with expensive (time consuming) objective functions, like our finite element analysis. The final form of this will be described elsewhere, for now we will show a hint of what is coming.

A possible approach is to fit a simple, easy to evaluate function to a modest number of objective function evaluations, and find the local minimum of this fitted function, rather than of the objective function itself. The fitted (or interpolated) function is subsequently improved in accuracy, by including data from the expensive function sampled only at the local minima of the fitted function. An improved estimate of the minimum is obtained on the refined interpolating function, and this cycle is repeated until a termination test is passed. This approach can be implemented in numerous ways, here we show our first attempt. The method chosen for MathOptimizer's next release uses a different fitting procedure, and it can handle any number of unknowns. This simple prototype can only handle one. The function optimize shown below uses Mathematica's Interpolation routine to construct a locally fitting parabolic interpolant, and the minimum of this is found with FindMinimum. Note that this approach gives a good solution to the tuning fork problem with just 16 function evaluations.

Clear[optimize];
optimize::out="search exceeded bounds";
Options[optimize]={MaxIterations->100};
optimize[obj_,{a_,b_},opts___?OptionQ]:=Catch[
Module[{opt,data,est,tol=0.05,result,prev={0,x->Null},maxit},(
maxit=MaxIterations/.Flatten[{opts}]/.Options[optimize];
data=N[Map[{#,obj[#]}&,Range[a,b,(b-a)/3]]];
opt=Interpolation[data,InterpolationOrder->2];   
result= FindMinimum[opt[x],{x,(a+b)/2,a,b},opts];
est=x/.result[[2]];
If[a > est || b < est,Throw[Message[optimize::out];prev],Null];
data=Union[Sort[Append[data,{est,obj[est]}]]];
While[Abs[prev[[1]]-result[[1]]] > tol && Length[data] < maxit,(
prev=result;
opt=Interpolation[data,InterpolationOrder->2];   
result= FindMinimum[opt[x],{x,est,a,b},opts];
est=x/.result[[2]];
If[a > est || b < est,Throw[Message[optimize::out];prev],Null];
data=Union[Sort[Append[data,{est,obj[est]}]]];
)];
result)]];

evals ={};
Timing[optimize[obj,{0.05,.15},AccuracyGoal->2,PrecisionGoal->2] ] 
Length[evals]

Do We Need Finite Elements For Such A Simple Shape?

One might guess that the lateral vibrations of tuning forks are simple enough to be studied analytically: [8] quotes an expression for the lowest natural frequency (Hz) of a tuning fork of height h and width of tines w (meters), YoungsModulus (Pascals) and Density (kg/(m^3)) as:

Clear[h, w, YoungsModulus, Density];
(Subscript[f, 1]) = ((π w (1.194^2) (Sqrt[(YoungsModulus/Density)]))/(16 (Sqrt[3]) (h^2)));

With such a simple relation no optimization would be needed, we could just solve this expression for h (here taking the positive root).

Solve[((Subscript[f, 1]) /. {YoungsModulus -> 190.*^9, Density -> 7700. , 
                w ->  .009}) == 440.0, h][[2]]

We see the height returned by the analytic formula is 35 % higher than our finite element results! This is because the analytic expression assumes that the tuning fork is built using bars terminated in a perfectly fixed end condition, whereas in our tuning fork the tines terminate in circular arcs, which are more compliant, leading to a lower frequency of resonance.

In case there is any doubt that the Mavart3D finite element result is accurate, here is the finite element calculation of the lowest natural frequency of a fixed-free beam that corresponds to the analytic expression above, and it gives 442.1 Hz, in reasonable agreement with the 440 Hz expected for a height of 0.128144 m. This could be improved slightly by refining the mesh, but perfect agreement is not possible because the analytic expression is for a thin beam, and so it does not contain Poisson's ratio, which is included in the finite element calculation.

\(\(\(\(h\  = 0.128144;\)

    \(w =  .009;\)

    Point[Delete]; \  (*\ Dimensions\ in\ meters\ *) 

    Point[New, {{0, 0}, {0, w}}];\ 

    \(Curve[New, {2, 1}];\)

    \(Curve[Translate, {1}, {w, 0}, 1];\)

    \(Surface[New, {1, 2}];\)

    \(Surface[Translate, {1}, {0, 0, h\ }, 1];\)

    \(Solid[New, {1, 2}];\)

    Node[Delete]; Mesh[2, 1, {2, 2, 10}];\ \ 

    \ \(Material[Delete];\)

    \(Material[New, 3, {7700. , 0, 0, 0, 190.*^9, 0.3}, "Steel"];\)\ \ 

    \(Fixity[Delete];\)

    \(Fixity[New, 1, 0.0, Node[Near, {\(,\)\(,\)\(0\)}]];\)

    \(Fixity[New, 2, 0.0, Node[Near, {\(,\)\(,\)\(0\)}]];\)

    \(Fixity[New, 3, 0.0, Node[Near, {\(,\)\(,\)\(0\)}]];\)

    \ \(WriteModel["mav.dat", 
        TITLE -> 
          StringJoin["Fixed-Free bar ", 
            ToString[Length[NodeList]], " Nodes "], 
        MACROS -> {"PRNT EIGE"}, RUNTASK -> {"EIGE"}, 
        NF -> 1];\)

    \(\(Analyze[]\) \); \  (*\ 
      Run\ the\ analysis\ and\ then\ retrieve\ the\ first\ eigen\ frequency\ \((Hz)\)\ *) 
\ \(ReadResult["mav.out", 
        ResultType -> "EigenValues"]\)\ [\([1, 
        2]\)]\)\(\ \)\)\)

This example shows that the classical textbook formula derived by Rayleigh [9] could not be used to tune the fork, and thus even for this apparently simple problem finite element analysis is required.

Conclusions

This paper described two Mathematica packages: MathOptimizer, a general purpose nonlinear (global and local) optimizer, and ModelMaker, a parametric finite element modeller. As an example, we presented the optimized design of a tuning fork, where the design objective function involves a finite element calculation. This simple example illustrates a typical engineering optimization problem where the objective function may be a black box (containing an elaborate calculation that is costly to evaluate). A prototype optimizer for expensive objective functions was demonstrated; an improved implementation of this will be included in the next release of MathOptimizer (and perhaps also offered as a separate application package).

Acknowledgements

This work was funded in part by Defence R&D Canada under Contract W7707-01-0746/001/HAL.

References

1. Pintér, J.D. Global Optimization in Action. Kluwer Academic Publishers, Dordrecht, (1996). 
2. Pintér, J.D.  Computational Global Optimization in Nonlinear Systems: An Interactive Tutorial. Lionheart Publishing, Inc., Atlanta, GA, (2001).
3. Pintér, J.D. MathOptimizer—An Advanced Modeling and Optimization System for Mathematica Users. User Guide. Pintér Consulting Services, Inc., Halifax, NS, (2002).
4. Trefethen, N.L. (2002) A Hundred-dollar, Hundred-digit Challenge. SIAM News, January-February 2000 issue, p. 3. See also for example, http://www.siam.org/siamnews/01-02/challenge.pdf and http://mathworld.wolfram.com/Hundred-DollarHundred-DigitChallengeProblems.html
5. Purcell, C.J., \(\(Building\ Finite\ Element\ Models\ with\ Mathematica\)\(,\)\) Mathematica Developers Conference, (1997).
6. Purcell, C.J., Modelling, 
          Analysis and Prototyping for Rapid Manufacturing,  Mathematica Developers Conference, (2001).
7. Wozniak, Z. and Purcell, C., MAVART 3D - A Review of Recent Developments, Proceedings of CanSmart, (2000). 
8. Bates, L.,  Beach. T.,  and Arnott, M., Determination of the Temperature Dependence of Young's Modulus for Stainless Steel Using a Tuning Fork, The Journal of Undergraduate Research in Physics, Vol. 18, No.1 (2002).
9. J.W.S. Rayleigh, The Theory of Sound, Dover (1878).