Distributed Solution of Optimal Control Problems Governed by Parabolic Equations
Matthias Heinkenschloss and Michael Herty

M. Heinkenschloss
Department of Computational and Applied Mathematics
Rice University

M. Herty
Fachbereich Mathematik, Technische Universität Kaiserslautern

We present a spatial domain decomposition (DD) preconditioner for the solution of discretized parabolic linear--quadratic optimal control problems. Our DD preconditioners are extensions of Neumann-Neumann DD preconditioners, which have been successfully applied to the solution of single elliptic partial differential equations and of linear--quadratic optimal control problems governed by elliptic equations.

The DD preconditioner is based on a decomposition of the spatial domain into non-overlapping subdomains. The optimality conditions for the parabolic linear quadratic optimal control problem is split into smaller problems restricted to spatial subdomain-time cylinders. These subproblems correspond to parabolic linear--quadratic optimal control problems on subdomains with Dirichlet data on interfaces. The coupling of these subdomain problems leads to a Schur complement system in which the unknowns are the state and adjoint variables on the subdomain interfaces in space and time.

The Schur complement system is solved using a preconditioned GMRES. The preconditioner is obtained from the solution of appropriate subdomain parabolic linear--quadratic optimal control problems. The dependence of the performance of these preconditioners on mesh size and subdomain size is studied numerically. Our tests indicate that their dependence on mesh size and subdomain size is similar to that of its counterpart applied to elliptic equations only. Our tests also suggest that the preconditioners are insensitive to the size of the control regularization parameter.

Key words: Optimal control, parabolic equations, domain decomposition, Neumann-Neumann methods.

December 2004

Preliminary Report on the Design and Implementation of Adifor90
Mike Fagan

In order to accurately and efficiently compute derivatives, many scientists and are abandoning divided differences in favor of Automatic Differentiation (AD). This paper describes the design and implementation of Adifor90, an AD tool for differentiating Fortran 90 programs. This paper also gives preliminary results that come from applying Adifor90 to UbikSolve, a component of the Truchas system.

Key words:

November 2004

Association-by-name and Association-by-address for Modern Programming Languages
Mike Fagan

Automatic Differentiation (AD) tools can employ 2 different techniques for associating derivatives with program variables. These 2 techniques are called association-by-address and association-by-name. Association-by-address is typically used for C-like languages. Association-by-name is typically used for Fortran-like languages. This report shows how to do association-by-name for C-like languages and association-by-address for Fortran-like languages.

Key words:

November 2004

Preliminary Report on the Design and Implementation of Adifor90
Mike Fagan

In order to accurately and efficiently compute derivatives, many scientists and are abandoning divided differences in favor of Automatic Differentiation (AD). This paper describes the design and implementation of Adifor90, an AD tool for differentiating Fortran 90 programs. This paper also gives preliminary results that come from applying Adifor90 to UbikSolve, a component of the Truchas system.

Key words:

November 2004

Activity Analysis in Adifor: Algorithms and Effectiveness
Mike Fagan and Alan Carle

Performance of automatic differentiation-enhaned codes may be improved by a source code analysis technique called activity analysis. Activity analysis seeks to limit the generated derivative code to the set of variables of interest to the programmer. This report describes both the static analysis and the runtime analysis methods.

Key words:

November 2004

Domain Decomposition Preconditioners for Linear-Quadratic Elliptic Optimal Control Problems
Matthias Heinkenschloss and Hoang Nguyen

We develop and analyze a class of overlapping domain decomposition (DD) preconditioners for linear-quadratic elliptic optimal control problems. Our preconditioners utilize the structure of the optimal control problems. Their execution requires the parallel solution of subdomain linear-quadratic elliptic optimal control problems, which are essentially smaller subdomain copies of the original problem. This work extends to optimal control problems the application and analysis of overlapping DD preconditioners, which have been used successfully for the solution of single PDEs.

We prove that the performance of the two-level versions of our preconditioners is independent of the mesh size and of the subdomain size. Our numerical studies indicate that the performance of our preconditioners for optimal control problems is comparable to the performance of their counterparts applied to single PDEs. Moreover, the performance of our preconditioners seems to be rather insensitive to the size of the control regularization parameter.

Key words: Domain decomposition, optimal control, preconditioner, KKT systems.

November 2004

The Standard Vector Library: A Software Framework for Coupling Complex Simulation and Optimizations
Anthony Padula, Shannon D. Scott, and William W. Symes

Object oriented design solves a fundamental programming problem arising in simulation driven optimization: the separation in code of multiple levels of abstraction naturally appearing in solution algorithms for such problems. The Standard Vector Library provides classes expressing core concepts (vector, function,...) of calculus in Hilbert space with minimal implementation dependence, and standardized interfaces behind which to hide application-dependent implementation details (data containers, function objects). Important innovations introduced by this project and its predecessor (the Hilbert Class Library) include vector space and function evaluation objects, natural product structures, and extensive tool or wrapper classes to ease application construction. The library features extensive use of ISO standard C++ support for both object-oriented and generic programming models, a component-friendly structure for support of distributed computing via client-server frameworks, and systematic extensibility of class capabilities through method-forwarding.

Key words:

September 2004

A New Algorithm for Continuation and Bifurcation Analysis of Large Scale Free Surface Flows
Zenaida Castillo

This thesis presents a new algorithm to find and follow particular solutions of parameterized nonlinear systems. Important applications often arise after spatial discretization of time dependent PDEs. We embed a block eigenvalue solver in a continuation framework for the computation of some specific eigenvalues of large Jacobian matrices that depend on one or more parameters. The new approach is then employed to study the behavior of an industrial process referred to as coating. Stability analysis of the discretized system that models this process is important because it provides alternatives for changing parameters in order to improve the quality of the final product or to increase productivity.

Experiments on several problems show the reliability of the new approach in the accurate detection of critical points. Further analysis of two-dimensional coating flow problems reveals that computational results are competitive with those of previous continuation approaches. As a byproduct, one obtains information about the stability of the process with no additional cost. Due to the size and structure of the matrices generated in three- dimensional free surface flow applications, it is necessary to use a general iterative linear solver, such as GMRES. However, GMRES displays a very slow rate of convergence as a consequence of the poor conditioning in the coefficient matrices. To speed up GMRES convergence, we developed and implemented a scalable approximate sparse inverse preconditioner. Numerical experiments demonstrate that this preconditioner greatly improves the convergence of the method. Results illustrate the effectiveness of the preconditioner on very large free surface flow problems with more than a million unknowns.

Key words:

August 2004 (Submitted in September 2004)

Gramians of Structured Systems and an Error Bound for Structure-Preserving Model Reduction
D.C. Sorensen & A.C Antoulas

In this paper a general framework is posed for defining the reachability and controllability gramians of structured linear dynamical systems. The novelty is that a formula for the gramian is given in the frequency domain. This formulation is surprisingly versatile and may be applied in a variety of structured problems. Moreover, this formulation enables a rather straightforward development of apriori error bounds for model reduction in the $\cH_2$ norm. The bound applies to a reduced model derived from projection onto the dominant eigenspace of the appropriate gramian. The reduced models are structure preserving because they arise as direct reduction of the original system in the reduced basis. A derivation of the bound is presented and verified computationally on a second order system arising from structural analysis.

Key words:

September 2004

Domain Decomposition Methods for Linear-Quadratic Elliptic Optimal Control Problems
Hoang Q. Nguyen

This thesis is concerned with the development of domain decomposition (DD) based preconditioners for linear-quadratic elliptic optimal control problems (LQ-EOCPs), their analysis, and numerical studies of their performance on model problems. The solution of LQ-EOCPs arises in many applications, either directly or as subproblems in Newton or Sequential Quadratic Programming methods for the solution of nonlinear elliptic optimal control problems. After a finite element discretization, convex LQ-EOCPs lead to large scale symmetric indefinite linear systems. The solution of these large systems is a very time consuming step and must be done iteratively, typically with a preconditioned Krylov subspace method. Developing good preconditioners for these linear systems is an important part of improving the overall performance of the solution method.

The DD preconditioners for LQ-EOCPs studied in this thesis are extensions of overlapping and nonoverlapping Neumann-Neumann DD preconditioners applied to single elliptic partial differential equations (PDEs). In our case, DD is applied on the optimization level. In particular, the proposed preconditioners require the parallel solution of subdomain optimal control problems that are related to restrictions of the original LQ-EOCP to a subdomain. Numerical results on several test problems have shown that the new preconditioners are effective. Their performance relative to decreases in finite element mesh size or increase in number of subdomains seem to be numerically comparable to that of overlapping and Neumann-Neumann preconditioners for single PDEs. Remarkably, the proposed preconditioners seem to be rather insensitive to control regularization parameters. For overlapping methods, theoretical results are provided to support the numerical observations.

Key words:

August 2004

An Infeasible Primal-Dual Algorithm for TV-Based Inf-Convolution-Type Image Restoration
M. Hintermueller & G. Stadler

In this paper, a primal-dual algorithm for TV-type image restoration is analyzed and tested. Analytically it turns out that employing a global L^s-regularization, with s>1, in the dual problem results in a local smoothing of the TV-regularization term in the primal problem. The local smoothing can alternatively be obtained as the infimal convolution of the l_r-norm, with r^{-1}+s^{-1}=1, and a smooth function. In the case r=s=2, this results in Gauss-TV-type image restoration. The globalized primal-dual algorithm introduced in this paper works with generalized derivatives, converges locally at a superlinear rate and is stable with respect to noise in the data. In addition, it utilizes a projection technique which reduces the size of the linear system that has to be solved per iteration. A comprehensive numerical study ends the paper.

Key words: Fenchel duality, generalized Newton-type methods, image restoration, total bounded variation.

August 2004

Shape Optimization in Unsteady Blood Flow: A Numerical Study of Non-Newtonian Effects
Feby Abraham, Marek Behr, & Matthias Heinkenschloss

This paper presents a numerical study of non-Newtonian effects on the solution of shape optimization problems involving unsteady pulsatile blood flow. We consider an idealized two-dimensional arterial graft geometry. Our computations are based on the Navier-Stokes equations generalized to non-Newtonian fluid, with the Carreau-Yasuda model employed to account for the shear-thinning behavior of blood. Using a gradient-based optimization algorithm, we compare the optimal shapes obtained using both the Newtonian and generalized Newtonian constitutive equations. Depending on the shear rate prevalent in the domain, substantial differences in the flow as well as in the computed optimal shape are observed when the Newtonian constitutive equation is replaced by the Carreau-Yasuda model. By varying a geometric parameter in our test case, we investigate the influence of the shear rate on the solution.

Key words: Shape optimization, blood flow, non-Newtonian fluids .

August 2004

A General Robust-Optimization Formulation for Nonlinear Programming
Yin Zhang

Most research in robust optimization has so far been focused on inequality-only, convex conic programming with simple linear models for uncertain parameters. Many practical optimization problems, however, are nonlinear and non-convex. Even in linear programming, coefficients can still be nonlinear functions of uncertain parameters. In this paper, we propose a formulation to extend the robust-optimization approach to a general nonlinear programming setting.

Key words: Parameterized nonlinear program, Robust optimization formulation.

July 2004

A Geometric Approach to Fluence Map Optimization in IMRT Cancer Treatment Planning
Yin Zhang and Michael Merritt

Intensity-modulated radiation therapy (IMRT) is a state-of-the-art technique for administering radiation to cancer patients. The goal of a treatment is to deliver a prescribed amount of radiation to the tumor, while limiting the amount absorbed by the surrounding healthy and critical organs. Planning an IMRT treatment requires determining fluence maps, each consisting of hundreds or more beamlet intensities. Since it is difficult or impossible to deliver a sufficient dose to a tumor without irradiating nearby critical organs, radiation oncologists have developed guidelines to allow tradeoffs by introducing so-called dose-volume constraints (DVCs), which specify a given percentage of volume for each critical organ that can be sacrificed if necessary. Such constraints, however, are of combinatorial nature and pose significant challenges to the fluence map optimization problem.

In this paper, we describe a new geometric approach to the fluence map optimization problem. Contrary to the traditional view, we regard dose distributions as our primary independent variables, while treating beamlet intensities as secondary ones. We consider two sets in the dose space: (i) the physical set consisting of physically realizable dose distributions, and (ii) the prescription set consisting of dose distributions meeting the prescribed tumor doses and satisfying the given dose-volume constraints. We seek a suitable dose distribution by successively projecting between these two sets. A crucial observation is that the projection onto the prescription set, which is non-convex, can be properly defined and easily computed. The projection onto the physical set, on the other hand, requires solving a nonnegative least squares problem. We show that this alternating projection algorithm is actually equivalent to a greedy algorithm driven by local sensitivity information readily available in our formulation. Moreover, the availability of such local sensitivity information will enable us to devise greedy algorithms to search for a desirable plan even when a ``good and achievable'' prescription is unknown.

Key words: Cancer radiation therapy, Optimal treatment planning, Fluence map optimization, Geometric Approach.

July 2004

Path-Following Methods for a Class of Constrained Minimization Problems in Function Space
M. Hintermueller and K. Kunisch

Path-following methods for primal-dual active set strategies requiring a regularization parameter are introduced. Existence of a path and its differentiability properties are analyzed. Monotonicity and convexity of the primal-dual path value function are investigated. Both feasible and infeasible approximations are considered. Numerical path following strategies are developed and their efficiency is demonstrated by means of examples.

Key words: Semi-smooth Newton methods, path-following methods, active set strategy, primal-dual methods.

July 2004

A Combined Shape-Newton and Topology Optimization Technique in Real-Time Image Segmentation
M. Hintermueller

In this paper, for solving a class of shape optimization problems, a new algorithmic concept combining shape and topological sensitivities is presented. The geometry of interest is represented by means of geometrical implicit functions. While the topology optimization phase is based on a gradient descent scheme, the phase based on shape sensitivity uses Newton-type descent flows. The algorithm is applied to edge detector based image segmentation problems, but it can be extended to solve more general shape optimization problems as well.

Key words: Image segmentation, level set method, shape optimization, shape sensitivity, shape Newton method, topological sensitivity.

July 2004

Filter Pattern Search Algorithms for Mixed Variable Constrained Optimization Problems
Mark A. Abramson, Charles Audet, and J.E. Dennis, Jr.

A new class of algorithms for solving nonlinearly constrained mixed variable optimization problems is presented. This class combines and extends the Audet-Dennis Generalized Pattern Search (GPS) algorithms for bound constrained mixed variable optimization, and their GPS-filter algorithms for general nonlinear constraints. In generalizing existing algorithms, new theoretical convergence results are presented that reduce seamlessly to existing results for more specific classes of problems. While no local continuity or smoothness assumptions are required to apply the algorithm, a hierarchy of theoretical convergence results based on the Clarke calculus is given, in which local smoothness dictate what can be proved about certain limit points generated by the algorithm. To demonstrate the usefulness of the algorithm, the algorithm is applied to the design of a load-bearing thermal insulation system. We believe this is the first algorithm with provable convergence results to directly target this class of problems.

Key words: Optimization, pattern search algorithm, filter methods, mixed variable programming, categorical variables, nonsmooth optimization.

June 2004

An Interior-Point Gradient Method for Large-Scale Totally Nonnegative Least Squares Problems
Michael Merritt and Yin Zhang

We study an interior-point gradient method for solving a class of so-called totally nonnegative least squares problems. At each iteration, the method decreases the residual norm along a diagonally scaled negative gradient direction with a special scaling. We establish the global convergence of the method, and present some numerical examples to compare the proposed method with some existing methods including the affine scaling method.

Key words: Interior-point gradient method, totally nonnegative least squares problem, global convergence.

May 2004

A Study of University Timetabling that Blends Graph Coloring with the Satisfaction of Various Essential and Preferential Conditions
Timothy A. Redl

Constructing a satisfactory conflict-free semester-long timetable of courses and creating a similarly satisfactory conflict-free timetable for end-of-semester final examinations are two closely related and often difficult problems that colleges and universities face each semester.

We discuss the relevance of such timetabling problems as a natural and practical application of graph coloring, and develop a mathematical and computational model for solving university timetabling problems using techniques of graph coloring that incorporates the satisfaction of both "essential" timetabling conditions (i.e., conditions or constraints that must be satisfied in order to produce a legal or feasible timetable) as well as suggested "preferential" timetabling conditions (i.e., additional conditions or constraints that need not necessarily be satisfied to produce a legal or legitimate timetable, but if satisfied may very well produce a more ``acceptable'' timetable for students and/or faculty members). We discuss in detail the step-by-step process that is taken to implement our timetabling-by-graph-coloring procedure, from the assembling of university course data, to creating a course conflict graph based on the assembled data, to coloring the conflict graph, to transforming this coloring to a conflict-free timetable, to finally assigning courses to classrooms. Once a conflict-free timetable of courses has been constructed, we present ways in which such a course timetable for a particular semester can be used to construct a conflict-free timetable of final examinations. Our model also considers a number of sociological scheduling concerns and preferences addressed by university registrars, faculty, staff, and students. Computational results, obtained by the author using actual data provided by Rice University and the University of St. Thomas, are documented.

May 2004

Interior-Point Gradient Methods with Diagonal-Scalings for Simple-Bound Constrained Optimization
Yin Zhang

In this paper, we study diagonally scaled gradient methods for simple-bound constrained optimization in a framework almost identical to that for unconstrained optimization, except that iterates are kept within the interior of the feasible region. We establish a satisfactory global convergence theory for such interior-point gradient methods applied to Lipschitz continuously differentiable functions without any further assumption. Moreover, a strong convergence result is obtained for a class of so-called L-nonlinear functions introduced in this paper which includes virtually all nonlinear functions that do not contain linear pieces.

Key words: Interior point gradient method, consistent scaling, L-nonlinear function, global convergence.

April 2004

Optimal Aeroacoustic Shape Design Using the Surrogate Management Framework
Alison L. Marsden, Meng Wang, J.E. Dennis, and Parviz Moin

Shape optimization is applied to time-dependent trailing-edge flow in order to minimize aerodynamic noise. Optimization is performed using the surrogate management framework (SMF), a non-gradient based pattern search method chosen for its efficiency and rigorous convergence properties. Using SMF, design space exploration is performed not with the expensive actual function but with an inexpensive surrogate function. The use of a polling step in the SMF guarantees that the algorithm generates a convergent subsequence of mesh points, each iterate of which is a local minimizer of the cost function on a mesh in the parameter space. Results are presented for an unsteady laminar flow past an acoustically compact airfoil. Constraints on lift and drag are handled within SMF by applying the filter pattern search method of Audet and Dennis, within which a penalty function is used to form and optimize a surrogate function. Optimal shapes that minimize noise have been identified for the trailing-edge problem in constrained and unconstrained cases. Results show a significant reduction (as much as 80%) in acoustic power with reasonable computational cost using several shape parameters. Physical mechanisms for noise reduction are discussed.

Key words: Surrogate optimization, derivative-free optimization, aeroacoustics, trailing-edge noise, pattern search, fluid mechanics

April 2004

Shape Optimization in Stationary Blood Flow: A Numerical Study of Non-Newtonian Effects
Feby Abraham, Marek Behr, and Matthias Heinkenschloss

We investigate the influence of the fluid constitutive model on the outcome of shape optimization tasks, motivated by optimal design problems in biomedical engineering. Our computations are based on the Navier-Stokes equations generalized to non-Newtonian fluid, with the Carreau-Yasuda model employed to account for the shear-thinning behavior of blood. The generalized Newtonian treatment exhibits striking differences in the velocity field for smaller shear rates. We apply sensitivity-based optimization procedure to a benchmark problem of flow through a right-angle cannula, and to a flow through an idealized arterial graft. For each of these problems we study the influence of the inflow velocity, and thus the shear rate. Furthermore for the arterial graft problem, we introduce an additional factor in the form of a geometric parameter, and study its effect on the optimal shape obtained.

Key words: Shape optimization, blood flow, non-Newtonian fluids .

February 2004

Second Order Behavior of Pattern Search Algorithms
Mark A. Abramson

Previous analyses of pattern search algorithms for unconstrained and linearly constrained minimization have focused on proving convergence of a subsequence of iterates to a limit point satisfying either directional or first-order necessary conditions for optimality, depending on the smoothness of the objective function in a neighborhood of the limit point. Even though pattern search methods require no derivative information, we are able to prove some limited directional second-order results. Although not as strong as classical second-order necessary conditions, these results are stronger than the first order conditions that many gradient-based methods satisfy. Under fairly mild conditions, we can eliminate from consideration all strict local maximizers and an entire class of saddle points.

Key words: Nonlinear Programming, pattern search algorithm, derivative-free optimization, convergence analysis, second order optimality conditions.

January 2004

Mesh Adaptive Direct Search Algorithms for Constrained Optimization
Charles Audet and J.E. Dennis, Jr.

This paper introduces the Mesh Adaptive Direct Search (MADS) class of algorithms for nonlinear optimization. MADS extends the Generalized Pattern Search (GPS) class by allowing local exploration, called polling, in a dense set of directions in the space of optimization variables. This means that under certain hypotheses, including a weak constraint qualification due to Rockafellar, MADS can treat constraints by the extreme barrier approach of setting the objective to infinity for infeasible points and treating the problem as unconstrained. The main GPS convergence result is to identify limit points where the Clarke generalized derivatives are nonnegative in a finite set of directions, called refining directions. Although in the unconstrained case, nonnegative combinations of these directions spans the whole space, the fact that there can only be finitely many GPS refining directions limits rigorous justification of the barrier approach to finitely many constraints for GPS. The MADS class of algorithms extend this result; the set of refining directions may even be dense in Rⁿ, although we give an example where it is not.

We present an implementable instance of MADS, and we illustrate and compare it with GPS on some test problems. We also illustrate the limitation of our results with examples.

Key words: Mesh adaptive direct search algorithms (MADS), convergence analysis, constrained optimization, nonsmooth analysis, Clarke directives, hypertangent, contingent cone.

January 2004 (Updated in October 2004)

Neumann-Neumann Domain Decomposition Preconditioners for Linear-Quadratic Elliptic Optimal Control Problems
Matthias Heinkenschloss and Hoang Nguyen

We present a class of domain decomposition (DD) preconditioners for the solution of elliptic linear-quadratic optimal control problems. Our DD preconditioners are extensions of Neumann-Neumann DD preconditioners, which have been successfully applied to the solution of single partial differential equations.

The DD preconditioners are based on a decomposition of the optimality conditions for the elliptic linear quadratic optimal control problem into smaller subdomain optimality conditions with Dirichlet boundary conditions for the states and the adjoints on the subdomain interfaces. These subdomain optimality conditions are coupled through Neumann interface conditions for the states and the adjoints. This decomposition leads to a Schur complement system in which the unknowns are the state and adjoint variables on the subdomain interfaces. The Schur complement operator is the sum of subdomain Schur complement operators, the application of which is shown to correspond to the solution of subdomain elliptic linear quadratic optimal control problems, which are essentially smaller copies of the original optimal control problem. We show that, under suitable conditions, the application of the inverse of the subdomain Schur complement operators requires the solution of a subdomain elliptic linear quadratic optimal control problem with Neumann interface conditions for the state.

The subdomain Schur complement operators are analyzed in the variational setting of the problem as well as the algebraic setting obtained after a finite element discretization of the problem. Definiteness properties of the algebraic form of the (subdomain) Schur complement operator(s) are studied. Numerical tests show that the dependence of these preconditioners on mesh size and subdomain size is comparable to its counterpart applied to elliptic equations only. These tests also show that the preconditioners are insensitive to the size of the control regularization parameter.

Key words: Optimal control, preconditioners, domain decomposition, Neumann-Neumann methods.

January 2004

Go to 2003 reports