Noah G. Harding Chair and Professor, Department of Computational and Applied Mathematics
This paper extends the algorithm of Benner, Heinkenschloss, Saak, and Weichelt: An inexact low-rank Newton-ADI method for large-scale algebraic Riccati equations, Applied Numerical Mathematics, 2016, Vol. 108, pages 125–142, DOI: 10.1016/j.apnum.2016.05.006, to Riccati equations associated with Hessenberg index-2 Differential Algebratic Equation (DAE) systems. Such DAE systems arise, e.g., from semi-discretized, linearized (around steady state) Navier-Stokes equations. The solution of the associated Riccati equation is important, e.g., to compute feedback laws that stabilize the Navier-Stokes equations. Challenges in the numerical solution of the Riccati equation arise from the large-scale of the underlying systems, the algebraic constraint in the DAE system, and the fact that matrices arising in some subproblems may only be marginally stable. These challenges are met by a careful extension of the inexact low-rank Newton-ADI method to the case of DAE systems. A main ingredient in the extension to the DAE case is the projection onto the manifold of the algebraic constraints. In the algorithm, the equations are never explicitly projected, but the projection is only applied as needed. The performance of the algorithm is illustrated on a large-scale Riccati equation associated with the stabilization of Navier-Stokes flow around a cylinder.
@article {PBenner_MHeinkenschloss_JSaak_HKWeichelt_2018a, author = {P. Benner and M. Heinkenschloss and J. Saak and H. K. Weichelt}, fauthor = {Peter Benner and Matthias Heinkenschloss and Jens Saak and Heiko K. Weichelt}, title = {Efficient Solution of Large-Scale Algebraic {R}iccati Equations Associated with Index-2 {DAEs} via the Inexact Low-Rank {N}ewton-{ADI} Method}, JOURNAL = {arXiv:1804.01410v1}, YEAR = {2018}, MONTH = {April}, URL = {http://arxiv.org/abs/1804.01410v1}, }
This paper addresses issues that originate in the extension of the Loewner framework to compute reduced order models (ROMs) of so-called quadratic-bilinear systems. The latter arise in semi-discretizations of fluid flow problems, such as Burgers' equation or the Navier-Stokes equations. In the linear case, the Loewner framework is data-driven and constructs a ROM from measurements of the transfer function; it does not explicitly require access to the system matrices, which is attractive in many settings. Research on extending the Loewner framework to quadratic-bilinear systems is ongoing. This paper presents one extension and provides details of its implementation that allow application to large-scale problems. This extension is applied to Burgers' equation. Numerical results illustrate the potential of the Loewner framework, but also expose additional issues that need to be addressed to make it fully applicable. Possible approaches to deal with some of these issues are outlined.
@INCOLLECTION{ACAntoulas_IVGosea_MHeinkenschloss_2018a, author = {A. C. Antoulas and I. V. Gosea and M. Heinkenschloss}, fauthor = {Athanasios C. Antoulas and Ion Victor Gosea and Matthias Heinkenschloss}, title = {On the {L}oewner Framework for Model Reduction of {B}urgers' Equation}, editor = {R. King}, booktitle = {Active Flow and Combustion Control 2018}, publisher = {Springer-Verlag}, address = {Berlin, Heidelberg, New York}, pages = {}, year = 2018, note = {submitted} }
This paper introduces an efficient sequential application of reduced order models (ROMs) to solve linear-quadratic optimal control problems with initial value controls. The numerical solution of such a problem requires Hessian-times-vector multiplications, each of which requires solving a linearized state equation with initial value given by the vector and solving a second order adjoint equation. Projection based ROMs are applied to these differential equations to generate a Hessian approximation. However, in general no fixed ROM well-approximates the application of the Hessian to all possible vectors of initial data. To improve a basic ROM, Heinkenschloss and Jando: Reduced Order Modeling for Time-Dependent Optimization Problems with Initial Value Controls (2016) introduce an augmentation of the basic ROM by the right hand side of the optimality system. This augmented ROM substantially improves the accuracy of the computed control, but this accuracy may still not be enough. The proposed sequential application of the augmented ROM can compute an approximate control with the same accuracy as the one obtained using only the expensive full order model, but at a fraction of the cost.
@INCOLLECTION {MHeinkenschloss_DJando_2018b, author = {M. Heinkenschloss and D. Jando}, fauthor = {Matthias Heinkenschloss and D{\"o}rte Jando}, title = {Sequential Reduced Order Modeling for Time-Dependent Optimization Problems with Initial Value Controls}, editor = {W. Keiper and A. Milde and S. Volkwein}, booktitle = {Reduced-Order Modeling ({ROM}) for Simulation and Optimization}, pages = {73--98}, year = {2018}, publisher = {Springer Verlag}, doi = {10.1007/978-3-319-75319-5_4}, url = {https://doi.org/10.1007/978-3-319-75319-5_4} }
This paper presents a new reduced order model (ROM) Hessian approximation for linear-quadratic optimal control problems where the optimal control is the initial value. Such problems arise in parameter identification, where the parameters to be identified appear in the initial data, and also as subproblems in multiple shooting formulations of more general optimal control problems. The new ROM Hessians can provide a substantially better approximation than the underlying basic ROM approximation, and thus can substantially reduce the computing time needed to solve these optimal control problems. The computation of a Hessian vector product requires the solution of the linearized state equation with initial value given by the vector to which the Hessian is applied to, followed by the solution of the second order adjoint equation. Projection based ROMs of these two linear differential equations are used to generate the Hessian approximation. The challenge is that in general no fixed ROM well-approximates the application of the Hessian to all possible vectors of initial data. The new approach, after having selected a basic ROM, augments this basic ROM by one vector. This vector is either the right hand side or the vector of initial data to which the Hessian is applied to. Although the size of the ROM increases only by one, this new augmented ROM produces substantially better approximations than the basic ROM. The use of these ROM Hessians in a conjugate gradient (CG) method is analyzed.
@article {MHeinkenschloss_DJando_2018a, AUTHOR = {M. Heinkenschloss and D. Jando}, FAUTHOR = {Matthias Heinkenschloss and D{\"o}rte Jando}, TITLE = {Reduced Order Modeling for Time-Dependent Optimization Problems with Initial Value Controls}, JOURNAL = {SIAM J. Sci. Comput.}, FJOURNAL = {SIAM Journal on Scientific Computing}, VOLUME = {40}, YEAR = {2018}, NUMBER = {1}, PAGES = {A22-A51}, DOI = {10.1137/16M1109084}, URL = {https://doi.org/10.1137/16M1109084} }
This paper proposes and analyzes two reduced order model (ROM) based approaches
for the efficient and accurate evaluation of the Conditional-Value-at-Risk
(CVaR) of quantities of interest (QoI) in engineering systems with uncertain
parameters. CVaR is used to model objective or constraint functions in risk averse
engineering design and optimization applications under uncertainty.
Evaluating the CVaR of the QoI requires sampling in the tail of the QoI distribution and typically requires
many solutions of an expensive full order model (FOM) of the engineering system.
Our ROM approaches substantially reduce this computational expense.
Both ROM based approaches use Monte-Carlo (MC) sampling. The first approach replaces the
computationally expensive FOM by inexpensive ROMs.
The resulting CVaR estimation error is proportional to the ROM error in the so-called risk-region,
a small region in the space of uncertain system inputs.
The second approach uses importance sampling (IS). ROM samples are used
to estimate the risk region and to construct a biasing distribution. Few FOM samples are then
drawn from this biasing distribution. Asymptotically as the ROM error goes to zero, the ROM-IS
sampling reduces the variance by a factor 1-β << 1, where
β in (0,1) is the quantile level at which CVaR is computed.
Numerical experiments on a system of semilinear convection-diffusion-reaction equations
illustrate the performance of the approaches.
@TECHREPORT{MHeinkenschloss_BKramer_TTakhtaganov_KWillcox_2017a, author = {M. Heinkenschloss and B. Kramer and T. Takhtaganov and K. Willcox}, fauthor = {Matthias Heinkenschloss and Boris Kramer and Timur Takhtaganov and Karen Willcox}, title = {Conditional-Value-at-Risk Estimation via Reduced-Order Models}, institution = {Department of Computational and Applied Mathematics}, address = {Rice University}, year = 2017, }
This paper provides a rigorous framework for the numerical solution of shape optimization problems in shell structure acoustics using a reference domain framework. The structure is modeled with Naghdi shell equations, fully coupled to boundary integral equations on a minimally regular surface, permitting the formulation of three-dimensional radiation and scattering problems on a two-dimensional set of reference coordinates. Well-posedness of this model and Fréchet differentiability of the state with respect to the surface shape are proven. For a class of shape optimization problems existence of optimal solutions under slightly stronger surface regularity assumptions is estiablished. Finally, adjoint equations are used to efficiently compute derivatives of the radiated field with respect to large numbers of shape parameters, which allows consideration of a rich space of shapes, and thus, of a broad range of design problems. A numerical example is presented to illustrate the applicability of the theoretical results.
@article{HAntil_SHardesty_MHeinkenschloss_2017a, fauthor = {Harbir Antil and Sean Hardesty and Matthias Heinkenschloss}, author = {H. Antil and S. Hardesty and M. Heinkenschloss}, title = {Shape Optimization of Shell Structure Acoustics}, journal = {SIAM Journal on Control and Optimization}, volume = {55}, number = {3}, pages = {1347-1376}, year = {2017}, doi = {10.1137/16M1070633}, URL = {http://dx.doi.org/10.1137/16M1070633} }
This supplement contains details that had to be omitted from the main paper H. Antil, S. Hardesty and M. Heinkenschloss: Shape Optimization of Shell Structure Acoustics, SIAM Journal on Control and Optimization, 2017, Vol. 55 (3), pp. 1347-1376, because of page limitations.
@TECHREPORT{HAntil_SHardesty_MHeinkenschloss_2017b, author = {H. Antil and S. Hardesty and M. Heinkenschloss}, title = {Supplementary Materials: Shape Optimization of Shell Structure Acoustics}, institution = {Department of Computational and Applied Mathematics}, address = {Rice University}, year = 2016, note = {Availabe electronically at http://www.caam.rice.edu/$\sim$heinken/mh\_publications.html} }
This paper introduces and analyzes a new parallel-in-time gradient type method for the solution of convex linear-quadratic discrete-time optimal control (DTOC) problems. Each iteration of the classical gradient method requires the solution of the forward-in-time state equation followed by the solution of the backward-in-time adjoint equation to compute the gradient. To introduce parallelism, the time steps are split into N groups corresponding to time subintervals. At the time subinterval boundaries state and adjoint information from the previous iteration is used. On each time subinterval the forward-in-time state equation is solved, the backward-in-time adjoint equation is solved, gradient-type information is generated, and the control are updated. These computations can be performed in parallel across time subintervals. State and adjoint information at time subinterval boundaries is then exchanged with neighboring subintervals and the process is repeated. The resulting iteration can be interpreted as a so-called (2N-1)-part iteration scheme. Convergence of the new parallel-in-time gradient type method is proven for suitable step-sizes by showing that an associated block companion matrix has spectral radius less than one. The performance of the new method is demonstrated on a DTOC problem obtained from a discretization of a 3D parabolic optimal control problem. In this example nearly perfect speed-up is observed for moderate number of time subdomains. This speed-up due to time decomposition multiplies existing speed-up due to parallelization in the solution of state and adjoint equations.
@TECHREPORT{XDeng_MHeinkenschloss_2016a, author = {X. Deng and M. Heinkenschloss}, fauthor = {Xiaodi Deng and Matthias Heinkenschloss}, title = {A Parallel-in-Time Gradient-Type Method for Discrete Time Optimal Control Problems}, institution = {Department of Computational and Applied Mathematics}, address = {Rice University}, type = {Preprint}, note = {Available from \url{http://www.caam.rice.edu/$\sim$heinken}}, }
This paper improves the inexact Kleinman-Newton method for solving algebraic Riccati equations by incorporating a line search and by systematically integrating the low-rank structure resulting from ADI methods for the approximate solution of the Lyapunov equation that needs to be solved to compute the Kleinman-Newton step. A convergence result is presented that tailors the convergence proof for general inexact Newton methods to the structure of Riccati equations and avoids positive semi-definiteness assumptions on the Lyapunov equation residual, which in general do not hold for low-rank approaches. In the convergence proof of this paper, the line-search is needed to ensure that the Riccati residuals decrease monotonically in norm. In the numerical experiments, the line search can lead to substantial reduction in overall number of ADI iterations and, therefore overall computational cost.
@article {PBenner_MHeinkenschloss_JSaak_HKWeichelt_2016a, AUTHOR = {P. Benner and M. Heinkenschloss and J. Saak and H. K. Weichelt}, FAUTHOR = {Benner, Peter and Heinkenschloss, Matthias and Saak, Jens and Weichelt, Heiko K.}, TITLE = {An inexact low-rank {N}ewton-{ADI} method for large-scale algebraic {R}iccati equations}, JOURNAL = {Appl. Numer. Math.}, FJOURNAL = {Applied Numerical Mathematics. An IMACS Journal}, VOLUME = {108}, YEAR = {2016}, PAGES = {125–142}, DOI = {10.1016/j.apnum.2016.05.006}, URL = {http://dx.doi.org/10.1016/j.apnum.2016.05.006} }
This paper improves the trust-region algorithm with adaptive sparse grids introduced in our previous paper for the solution of optimization problems governed by partial differential equations (PDEs) with uncertain coefficients. The previous algorithm used adaptive sparse grid discretizations to generate models that are applied in a trust-region framework to generate a trial step. The decision whether to accept this trial step as the new iterate, however, required relatively high fidelity adaptive discretizations of the objective function. In this paper, we extend the algorithm and convergence theory to allow the use of low-fidelity adaptive sparse-grid models in objective function evaluations. This is accomplished by extending conditions on inexact function evaluations used in previous trust-region frameworks. Our algorithm adaptively builds two separate sparse grids: one to generate optimization models for the step computation, and one to approximate the objective function. These adapted sparse grids typically contain significantly fewer points than the high-fidelity grids, which leads to a dramatic reduction in the computational cost. This is demonstrated numerically using two examples. Moreover, the numerical results indicate that the new algorithm rapidly identifies the stochastic variables that are relevant to obtaining an accurate optimal solution. When the number of such variables is independent of the dimension of the stochastic space, the algorithm exhibits near dimension-independent behavior.
@article{DPKouri_MHeinkenschloss_DRidzal_BGvanBloemenWaanders_2014a, author = {D. P. Kouri and M. Heinkenschloss and D. Ridzal and B. G. {van Bloemen Waanders} }, title = {Inexact Objective Function Evaluations in a Trust-Region Algorithm for {PDE}-Constrained Optimization under Uncertainty}, journal = {SIAM Journal on Scientific Computing}, volume = {36}, number = {6}, pages = {A3011-A3029}, year = {2014}, doi = {10.1137/140955665}, url = {http://dx.doi.org/10.1137/140955665} }
We develop and analyze a trust-region sequential quadratic programming (SQP) method for the solution of smooth equality constrained optimization problems, which allows the inexact and hence iterative solution of linear systems. Iterative solution of linear systems is important in large-scale applications, such as optimization problems with partial differential equation constraints, where direct solves are either too expensive or not applicable. Our trust-region SQP algorithm is based on a composite-step approach that decouples the step into a quasi-normal and a tangential step. The algorithm includes critical modifications of substep computations needed to cope with the inexact solution of linear systems. The global convergence of our algorithm is guaranteed under rather general conditions on the substeps. We propose algorithms to compute the substeps and prove that these algorithms satisfy global convergence conditions. All components of the resulting algorithm are specified in such a way that they can be directly implemented. Numerical results indicate that our algorithm converges even for very coarse linear system solves.
@article{MHeinkenschloss_DRidzal_2014a, author = {M. Heinkenschloss and D. Ridzal}, title = {A Matrix-Free Trust-Region {SQP} Method for Equality Constrained Optimization}, journal = {SIAM Journal on Optimization}, volume = {24}, number = {3}, pages = {1507-1541}, year = {2014}, doi = {10.1137/130921738}, url = {http://dx.doi.org/10.1137/130921738} }
Projection based methods lead to reduced order models (ROMs) with dramatically reduced numbers of equations and unknowns. However, for nonlinear or parametrically varying problems the cost of evaluating these ROMs still depends on the size of the full order model and therefore is still expensive. The Discrete Empirical Interpolation Method (DEIM) further approximates the nonlinearity in the projection based ROM. The resulting DEIM ROM nonlinearity depends only on a few components of the original nonlinearity. If each component of the original nonlinearity depends only on a few components of the argument, the resulting DEIM ROM can be evaluated efficiently at a cost that is independent of the size of the original problem. For systems obtained from finite difference approximations, the i th component of the original nonlinearity often depends only on the ith component of the argument. This is different for systems obtained using finite element methods, where the dependence is determined by the mesh and by the polynomial degree of the finite element subspaces. This paper describes two approaches of applying DEIM in the finite element context, one applied to the assembled and the other to the unassembled form of the nonlinearity. We carefully examine how the DEIM is applied in each case, and the substantial efficiency gains obtained by the DEIM. In addition, we demonstrate how to apply DEIM to obtain ROMs for a class of parameterized system that arises, e.g., in shape optimization. The evaluations of the DEIM ROMs are substantially faster than those of the standard projection based ROMs. Additional gains are obtained with the DEIM ROMs when one has to compute derivatives of the model with respect to the parameter.
@incollection {HAntil_MHeinkenschloss_DCSorensen_2014a, AUTHOR = {H. Antil and M. Heinkenschloss and D. C. Sorensen}, TITLE = {Application of the Discrete Empirical Interpolation Method to Reduced Order Modeling of Nonlinear and Parametric Systems}, EDITOR = {A. Quarteroni and G. Rozza}, BOOKTITLE = {Reduced Order Methods for Modeling and Computational Reduction}, SERIES = {MS\&A. Model. Simul. Appl.}, VOLUME = {9}, PAGES = {101-136}, PUBLISHER = {Springer Italia, Milan}, YEAR = {2014}, DOI = {10.1007/978-3-319-02090-7\_4} }
"Appendix": Proof of the Local Discretization Error Estimate for the Optimal Control Problem in the Presence of Interior Layers.