Publications

This work analyzes the sensitivities of the solution of a system of ordinary differential equations (ODEs) and a corresponding quantity of interest (QoI) to perturbations in a state-dependent component function that appears in the governing ODEs. This extends existing ODE sensitivity results, which consider the sensitivity of the ODE solution with respect to state-independent parameters. It is shown that with Caratheodory-type assumptions on the ODEs, the Implicit Function Theorem can be applied to establish continuous Frechet differentiability of the ODE solution with respect to the component function. These sensitivities are used to develop new estimates for the change in the ODE solution or QoI when the component function is perturbed. In applications, this new sensitivity-based bound on the ODE solution or QoI error is often much tighter than classical Gronwall-type error bounds. The sensitivity-based error bounds are applied to Zermelo's problem and to a trajectory simulation for a hypersonic vehicle.

@article {JRCangelosi_MHeinkenschloss_2024a,
     AUTHOR = {J. R. Cangelosi and M. Heinkenschloss},
    FAUTHOR = {Jonathan R. Cangelosi and Matthias Heinkenschloss},
      TITLE = {Sensitivity of {ODE} Solutions with Respect to Component Functions in the Dynamics},
    JOURNAL = {arXiv:2411.09655v1},
       YEAR = {2024},
        URL = {https://doi.org/10.48550/arXiv.2411.09655},
        DOI = {10.48550/arXiv.2411.09655}
}

This paper develops a line-search algorithm that uses objective function models with tunable accuracy to solve smooth optimization problems with convex constraints. The evaluation of objective function and its gradient is potentially computationally expensive, but it is assumed that one can construct effective, computationally inexpensive models. This paper specifies how these models can be used to generate new iterates. At each iteration, the model has to satisfy function error and relative gradient error tolerances determined by the algorithm based on its progress. Moreover, a bound for the model error is used to explore regions where the model is sufficiently accurate. The algorithm has the same first-order global convergence properties as standard line-search methods, but only uses the models and the model error bounds. The algorithm is applied to problems where the evaluation of the objective requires the solution of a large-scale system of nonlinear equations. The models are constructed from reduced order models of this system. Numerical results for partial differential equation constrained optimization problems show the benefits of the proposed algorithm.

@article {DSGrundvig_MHeinkenschloss_2024a,
    AUTHOR = {D. S. Grundvig and M. Heinkenschloss},
   FAUTHOR = {Dane S. Grundvig and Matthias Heinkenschloss},
     TITLE = {Line-Search Based Optimization using Function Approximations with Tunable Accuracy},
   JOURNAL = {Optim. Methods Softw.},
  FJOURNAL = {Optimization Methods and Software},
      YEAR = {2024},
       DOI = {10.1080/10556788.2024.2436192},
       URL = {https://doi.org/10.1080/10556788.2024.2436192},
      NOTE = {online first}
}

A new diagonalization technique for the parallel-in-time solution of linear-quadratic optimal control problems with time-invariant system matrices is introduced. The target problems are often derived from a semi-discretization of a Partial Differential Equation (PDE)-constrained optimization problem. The solution of large-scale time dependent optimal control problems is computationally challenging as the states, controls, and adjoints are coupled to each other throughout the whole time domain. This computational difficulty motivates the use of parallel-in-time methods. For time-periodic problems our diagonalization efficiently transforms the discretized optimality system into $n_t$ (=number of time steps) decoupled complex valued $2n_y \times 2n_y$ systems, where $n_y$ is the dimension of the state space. These systems resemble optimality systems corresponding to a steady-state version of the optimal control problem and they can be solved in parallel across the time steps, but are complex valued. For optimal control problems with initial value state equations a direct solution via diagonalization is not possible, but an efficient preconditioner can be constructed from the corresponding time periodic optimal control problem. The preconditioner can be efficiently applied parallel-in-time using the diagonalization technique. The observed number of preconditioned GMRES iterations is small and insensitive to the size of the problem discretization.

@article {MHeinkenschloss_NJKroeger_2024a,
    AUTHOR = {M. Heinkenschloss and N. J. Kroeger},
   FAUTHOR = {Matthias Heinkenschloss and Nathaniel J. Kroeger},
     TITLE = {A new diagonalization based method for parallel-in-time solution
              of linear-quadratic optimal control  problems},
   JOURNAL = {ESAIM Control Optim. Calc. Var.},
  FJOURNAL = {ESAIM. Control, Optimisation and Calculus of Variations},
      YEAR = {2024},
    VOLUME = {30},
    NUMBER = {62}, %Article Number 62
       DOI = {10.1051/cocv/2024051},
       URL = {https://doi.org/10.1051/cocv/2024051}
}

This paper integrates nonlinear-manifold reduced order models (NM-ROMs) with domain decomposition (DD). NM-ROMs approximate the full order model (FOM) state in a nonlinear-manifold by training a shallow, sparse autoencoder using FOM snapshot data. These NM-ROMs can be advantageous over linear-subspace ROMs (LS-ROMs) for problems with slowly decaying Kolmogorov $n$-width. However, the number of NM-ROM parameters that need to be trained scales with the size of the FOM. Moreover, for ``extreme-scale" problems, the storage of high-dimensional FOM snapshots alone can make ROM training expensive. To alleviate the training cost, this paper applies DD to the FOM, computes NM-ROMs on each subdomain, and couples them to obtain a global NM-ROM. This approach has several advantages: Subdomain NM-ROMs can be trained in parallel, involve fewer parameters to be trained than global NM-ROMs, require smaller subdomain FOM dimensional training data, and can be tailored to subdomain-specific features of the FOM. The shallow, sparse architecture of the autoencoder used in each subdomain NM-ROM allows application of hyper-reduction (HR), reducing the complexity caused by nonlinearity and yielding computational speedup of the NM-ROM. This paper provides the first application of NM-ROM (with HR) to a DD problem. In particular, this paper details an algebraic DD reformulation of the FOM, training a NM-ROM with HR for each subdomain, and a sequential quadratic programming (SQP) solver to evaluate the coupled global NM-ROM. Theoretical convergence results for the SQP method and {\it a priori} and {\it a posteriori} error estimates for the DD NM-ROM with HR are provided. The proposed DD NM-ROM with HR approach is numerically compared to a DD LS-ROM with HR on the 2D steady-state Burgers' equation, showing an order of magnitude improvement in accuracy of the proposed DD NM-ROM over the DD LS-ROM.

@article {ANDiaz_YChoi_MHeinkenschloss_2024a,
    AUTHOR = {A. N. Diaz and Y. Choi and M. Heinkenschloss},
   FAUTHOR = {Alejandro N. Diaz and Youngsoo Choi and Matthias Heinkenschloss},
     TITLE = {A fast and accurate domain-decomposition nonlinear manifold reduced order model},
   JOURNAL = {Comput. Methods Appl. Mech. Engrg.},
  FJOURNAL = {Computer Methods in Applied Mechanics and Engineering},
    VOLUME = {425},
      YEAR = {2024},
     PAGES = {116943},
       DOI = {10.1016/j.cma.2024.116943},
       URL = {https://doi.org/10.1016/j.cma.2024.116943}
}

This manuscript describes a methodology for simultaneous vehicle and trajectory optimization of a hypersonic glide vehicle. The co-design problem is formulated as an optimization problem with constraints including vehicle dynamics, path constraints (e.g., surface heating), and other constraints. The discretized optimization problem is solved simultaneously in the vehicle design parameters, the state variables, and the controls using an interior point method. Gaussian process (GP) surrogates, which are generated from sample candidate designs and flight conditions, are used to model vehicle aerodynamic performance and mass properties, as well as their first and second-order derivatives required by the optimizer. These GP surrogates and their derivatives are computationally inexpensive, making the all-at-once optimization approach for the co-design problem more tractable. To mitigate the effect of surrogate model errors on the solution of the optimal control problem, the GP models are refined using samples of the vehicle aerodynamic performance and mass properties at the solution of the co-design problem with the current surrogate. The resulting framework is applied to maximizing the range of a hypersonic glide vehicle with path and terminal constraints. Possible extensions of this methodology are also discussed, including the incorporation of more complex vehicle models such as multi-fidelity models, as well as adaptive surrogate modeling strategies to mitigate the effect of model errors on the solution of the optimal control problem.

@inproceedings{JRCangelosi_MHeinkenschloss_JTNeedels_JJAlonso_2024a,
   fauthor = {Jonathan R. Cangelosi and Matthias Heinkenschloss and Jacob T. Needels
               and Juan J. Alonso},
    author = {J. R. Cangelosi and M. Heinkenschloss and J. T. Needels
               and J. J. Alonso},
     title = {Simultaneous Design and Trajectory Optimization for Boosted Hypersonic
              Glide Vehicles},
 booktitle = {Paper AIAA 2024-0375.
             2024 AIAA Science and Technology Forum and Exposition (AIAA SciTech Forum)},
      year = {2024},
       doi = {10.2514/6.2024-0375},
       url = {https://doi.org/10.2514/6.2024-0375}
}

A nonlinear-manifold reduced order model (NM-ROM) is a great way of incorporating underlying physics principles into a neural network-based data-driven approach. We combine NM-ROMs with domain decomposition (DD) for efficient computation. NM-ROMs offer benefits over linear-subspace ROMs (LS-ROMs) but can be costly to train due to parameter scaling with the full-order model (FOM) size. To address this, we employ DD on the FOM, compute subdomain NM-ROMs, and then merge them into a global NM-ROM. This approach has multiple advantages: parallel training of subdomain NM-ROMs, fewer parameters than global NM-ROMs, and adaptability to subdomain-specific FOM features. Each subdomain NM-ROM uses a shallow, sparse autoencoder, enabling hyper-reduction (HR) for improved computational speed. In this paper, we detail an algebraic DD formulation for the FOM, train HR-equipped NM-ROMs for subdomains, and numerically compare them to DD LS-ROMs with HR. Results show a significant accuracy boost, on the order of magnitude, for the proposed DD NM-ROMs over DD LS-ROMs in solving the 2D steady-state Burgers' equation.

@inproceedings{ANDiaz_YChoi_MHeinkenschloss_2023a, 
    AUTHOR = {A. N. Diaz and Y. Choi and M. Heinkenschloss},
   FAUTHOR = {Alejandro N. Diaz and Youngsoo Choi and Matthias Heinkenschloss},
     TITLE = {Nonlinear-manifold reduced order models with domain decomposition},
 BOOKTITLE = {Proceedings of the NeurIPS 2023 Workshop: Machine Learning and the Physical Sciences
              as part of the 37th International Conference on
              Neural Information Processing Systems (NeurIPS 2023)},
      YEAR = {2023},
      NOTE = {Also available as arXiv:2312.00713},
       DOI = {10.48550/arXiv.2312.00713},
       URL = {https://doi.org/10.48550/arXiv.2312.00713}
}

This paper extends interpolatory model reduction to systems with (up to) quadratic-bilinear dynamics and quadratic-bilinear outputs. These systems are referred to as QB-QB systems and arise in a number of applications, including fluid dynamics, optimal control, and uncertainty quantification. In the interpolatory approach, the reduced order models (ROMs) are based on a Petrov-Galerkin projection and the projection matrices are constructed so that transfer function components of the ROM interpolate the corresponding transfer function components of the original system. To extend the approach to systems with QB outputs, this paper derives system transfer functions and sufficient conditions on the projection matrices that guarantee the aforementioned interpolation properties. Alternatively, if the system has linear dynamics and quadratic outputs, one can introduce auxiliary state variables to transform it into a system with QB dynamics and linear outputs to which known interpolatory model reduction can be applied. This transformation approach is compared with the proposed extension that directly treats quadratic outputs. The comparison shows that transformation hides problem structure. Numerical examples illustrate that keeping the original QB-QB structure leads to ROMs with better approximation properties.

@article{ANDiaz_IVGosea_MHeinkenschloss_ACAntoulas_2023a,
    AUTHOR = {A. N. Diaz and I. V. Gosea and M. Heinkenschloss and A. C. Antoulas},
   FAUTHOR = {Alejandro N. Diaz and Ion Victor Gosea and
               Matthias Heinkenschloss and Athanasios C. Antoulas},
     TITLE = {Interpolation-Based Model Reduction of Quadratic-Bilinear Dynamical 
              Systems with Quadratic-Bilinear Outputs},
   JOURNAL = {Advances in Computational Mathematics},
    VOLUME = {49},
    NUMBER = {95},
      YEAR = {2023},
       DOI = {10.1007/s10444-023-10096-2},
       URL = {https://doi.org/10.1007/s10444-023-10096-2}
} 

This paper introduces reduced order model (ROM) based Hessian approximations for use in inexact Newton methods for the solution of optimization problems implicitly constrained by a large-scale system, typically a discretization of a partial differential equation (PDE). The direct application of an inexact Newton method to this problem requires the solution of many PDEs per optimization iteration. To reduce the computational complexity, a ROM Hessian approximation is proposed. Since only the Hessian is approximated, but the original objective function and its gradient is used, the resulting inexact Newton method maintains the first-order global convergence property, under suitable assumptions. Thus even computationally inexpensive lower fidelity ROMs can be used, which is different from ROM approaches that replace the original optimization problem by a sequence of ROM optimization problem and typically need to accurately approximate function and gradient information of the original problem. In the proposed approach, the quality of the ROM Hessian approximation determines the rate of convergence, but not whether the method converges. The projection based ROM is constructed from state and adjoint snapshots, and is relatively inexpensive to compute. Numerical examples on semilinear parabolic optimal control problems demonstrate that the proposed approach can lead to substantial savings in terms of overall PDE solves required.

@INCOLLECTION{MHeinkenschloss_CMagruder_2022a,
   AUTHOR = {M. Heinkenschloss and C. Magruder},
   FAUTHOR = {Matthias Heinkenschloss and Caleb Magruder},
   title = {Reduced Order Model {H}essian Approximations in {N}ewton Methods for Optimal Control},
   editor = {C. A. Beattie and P. Benner and M. Embree and S. Gugercin and S. Lefteriu},
   booktitle = {Realization and Model Reduction of Dynamical Systems -
                A Festschrift in Honor of the 70th Birthday of {T}hanos {A}ntoulas},
   publisher = {Springer International Publishing},
   address = {Cham},
   pages = {335--351},
   year = 2022,
   URL = {https://doi.org/10.1007/978-3-030-95157-3_18},
   DOI = {10.1007/978-3-030-95157-3_18}
}

The Loewner framework is extended to compute reduced order models (ROMs) for systems governed by the incompressible Navier-Stokes (NS) equations. For quadratic ordinary differential equations (ODEs) it constructs a ROM directly from measurements of transfer function components derived from an expansion of the system's input-to-output map. Given measurements, no explicit access to the system is required to construct the ROM. To extend the Loewner framework, the NS equations are transformed into ODEs by projecting onto the subspace defined by the incompressibility condition. This projection is used theoretically, but avoided computationally. This paper presents the overall approach. Currently, transfer function measurements are obtained via computational simulations; obtaining them from experiments is an open issue. Numerical results show the potential of the Loewner framework, but also reveal possible lack of stability of the ROM. A possible approach, which currently requires access to the NS system, to deal with these instabilities is outlined.

@INCOLLECTION{ANDiaz_MHeinkenschloss_2022a,
   author = {A. N. Diaz and M. Heinkenschloss},
   fauthor = {Alejandro N. Diaz and Matthias Heinkenschloss},
   title = {Towards Data-Driven Model Reduction of the {N}avier-{S}tokes Equations using the {L}oewner Framework},
   editor = {R. King and D. Peitsch},
   booktitle = {Active Flow and Combustion Control 2021},
   publisher = {Springer International Publishing},
   address = {Cham},
   year = 2022,
   pages = {225--239},
   DOI = {10.1007/978-3-030-90727-3_14},
   URL = {https://doi.org/10.1007/978-3-030-90727-3_14}
}

This book highlights new developments in the wide and growing field of partial differential equations (PDE)-constrained optimization. Optimization problems where the dynamics evolve according to a system of PDEs arise in science, engineering, and economic applications and they can take the form of inverse problems, optimal control problems or optimal design problems. This book covers new theoretical, computational as well as implementation aspects for PDE-constrained optimization problems under uncertainty, in shape optimization, and in feedback control, and it illustrates the new developments on representative problems from a variety of applications.

@book{RHerzog_MHeinkenschloss_DKalise_GStadler_ETrelat_2022a,
      editor = {R. Herzog and M. Heinkenschloss and D. Kalise and G. Stadler and E. Tr\'elat},
     feditor = {Roland Herzog and Matthias Heinkenschloss and Dante Kalise and Georg Stadler and Emmanuel Tr\'elat},
       title = {Optimization and Control for Partial Differential Equations:
                Uncertainty Quantification, Open and Closed-Loop Control, and Shape Optimization},
   publisher = {De Gruyter},
      series = {Radon Series on Computational and Applied Mathematicsi, Vol.~29},
        year = {2022}
}

Bioethanol is an important biofuel with high potential to be a safety environmental alternative to fossil fuels. However, bioethanol production involves many distillation steps, which is one of the most energy consuming separation processes. Unconventional distillation techniques were proposed, aiming to reduce energy cost and CO2 emissions on bioethanol production. Parastillation and metastillation processes and novelty combinations of these techniques were explored. In all proposed configurations the liquid or vapor phase was divided inside the column. For the first time, parastillation columns with multiple condensers and columns with different distillation techniques combined into one unique column shell, here denominated as combined columns, were studied. Multiple condensers set, in parastillation columns, avoids the non-necessary mixed of the vapor phases. Combined columns are composed of parastillation trays in the rectifying section and of metastillation or conventional trays in the stripping section. The new configurations reduce total annual costs and CO2 emissions up to 35% and 42%, respectively, when compared with traditional distillation. Reduction in the column diameter by using metastillation instead of conventional distillation was possible without an increase in the column height, differently from previously results. Therefore, biofuel and neutral alcohol production can be even more economical and sustainable.

@article{LCKBiasi_FRMBatista_RJZemp_ALRRomano_MHeinkenschloss_AJAMeirelles_2021a,
    AUTHOR = {L. C. K. Biasi and  F. R. M. Batista and  R. J. Zemp
              and A. L. R. Romano and  M. Heinkenschloss and  A. J. A. Meirelles},
   FAUTHOR = {Lilian C. Kramer Biasi and  Fabio R. M. Batista and  Roger J. Zemp
              and Ana Laura Ramos Romano and Matthias Heinkenschloss and  Antonio J. A. Meirelles},
     title = {Parastillation and metastillation applied to bioethanol and neutral alcohol purification
             with energy savings},
   journal = {Chemical Engineering and Processing - Process Intensification},
    volume = {162},
     pages = {108334},
      year = {2021},
       doi = {10.1016/j.cep.2021.108334},
       url = {https://doi.org/10.1016/j.cep.2021.108334}
}

This paper introduces a modified version of the recent data-driven Loewner framework to compute reduced order models (ROMs) for a class of semi-explicit differential algebraic equation (DAE) systems, which include the semi-discretized linearized Navier-Stokes /Oseen equations. The modified version estimates the polynomial part of the original transfer function from data and incorporate this estimate into the Loewner ROM construction. Without this proposed modification the transfer function of the Loewner ROM is strictly proper, i.e., goes to zero as the magnitude of the frequency goes to infinity, and therefore may have a different behavior for large frequencies than the transfer function of the original system. The modification leads to a Loewner ROM with a transfer function that has a strictly proper and a polynomial part, just as the original model. This leads to better approximations for transfer functioncomponents in which the coefficients in the polynomial part are not too small. The construction of the improved Loewner ROM is described and the improvement is demonstrated on a large-scale system governed by the semi-discretized Oseen equations.

@INCOLLECTION{ACAntoulas_IVGosea_MHeinkenschloss_2020a,
   author = {A. C. Antoulas and I. V. Gosea and M. Heinkenschloss},
   fauthor = {Athanasios C. Antoulas and Ion Victor Gosea and Matthias Heinkenschloss},
   title = {Data-Driven Model Reduction for a Class of Semi-Explicit
            {DAE}s Using the {L}oewner Framework},
   editor = {T. Reis and S. Grundel and S. Sch{\"o}ps},
   booktitle = {{P}rogress in {D}ifferential-{A}lgebraic {E}quations {II}},
   publisher = {Springer International Publishing},
   address = {Cham},
   year = 2020,
   pages = {185--210},
   doi = {10.1007/978-3-030-53905-4_7},
   url = {https://doi.org/10.1007/978-3-030-53905-4_7}
}
	  

This paper presents a new unified model that allows the simulation of distillation columns with any integer number of phase divisions. Such columns have the potential for substantial savings in energy and capital costs. However, conventional simulators do not cover columns with phase divisions and previous simulations required a tailored algorithm for each type of column. The proposed model uses a unique set of MESH-equations for parastillation, metastillation, and conventional distillation. With small modifications in the stage indexing, the model allows the simulation of conventional, top and bottom-DWCs. This generalization was possible by introducing two new variables: the number of liquid (θ) and vapor (β) phase divisions. The positive impact of increasing the number of phase divisions was demonstrated in the bioethanol distillation, by analyzing parastillation and metastillation columns. These columns reduce the operational and total annual costs in 19% and 15%, when compared to conventional columns.

@article{LCKBiasi_MHeinkenschloss_FRMBatista_RJZemp_ALRRomano_AJAMeirelles_2020a,
    AUTHOR = {L. C. K. Biasi and  M. Heinkenschloss and  F. R. M. Batista and  R. J. Zemp
              and A. L. R. Romano and  A. J. A. Meirelles},
   FAUTHOR = {Lilian C. Kramer Biasi and  Matthias Heinkenschloss and  Fabio R. M. Batista and  Roger J. Zemp
              and Ana Laura Ramos Romano and  Antonio J. A. Meirelles},
    title = {A new unified model to simulate columns with multiple phase divisions
             and their impact on energy savings},
    journal = {Computers \& Chemical Engineering},
    volume = {140},
    pages = {106937},
    year = {2020},
    doi = {10.1016/j.compchemeng.2020.106937},
    url = {https://doi.org/10.1016/j.compchemeng.2020.106937}
}

This paper shows how to systematically and efficiently improve a reduced-order model (ROM) to obtain a better ROM-based estimate of the Conditional Value-at-Risk (CVaR) of a computationally expensive quantity of interest (QoI). Efficiency is gained by exploiting the structure of CVaR, which implies that a ROM used for CVaR estimation only needs to be accurate in a small region of the parameter space, called the epsilon risk region. Hence, any full-order model (FOM) queries needed to improve the ROM can be restricted to this small region of the parameter space, thereby substantially reducing the computational cost of ROM construction. However, an example is presented which shows that simply constructing a new ROM that has a smaller error with the FOM is in general not sufficient to yield a better CVaR estimate. Instead a combination of previous ROMs is proposed that achieves a guaranteed improvement, as well as epsilon risk regions that converge monotonically to the FOM risk region with decreasing ROM error. Error estimates for the ROM-based CVaR estimates are presented. The gains in efficiency obtained by improving a ROM only in the small epsilon risk region over a traditional greedy procedure on the entire parameter space is illustrated numerically.

@article {MHeinkenschloss_BKramer_TTakhtaganov_2020a,
    AUTHOR = {M. Heinkenschloss and B. Kramer and T. Takhtaganov},
   FAUTHOR = {Matthias Heinkenschloss and Boris Kramer and
              Timur Takhtaganov},
     TITLE = {Adaptive Reduced-Order Model Construction for Conditional Value-at-Risk Estimation
   JOURNAL = {SIAM/ASA J. Uncertainty Quantification},
  FJOURNAL = {SIAM/ASA Journal on Uncertainty Quantification},
    VOLUME = {8},
      YEAR = {2020},
    NUMBER = {2},
     PAGES = {668-692},
       DOI = {10.1137/19M1257433},
       URL = {https://doi.org/10.1137/19M1257433}
}

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 and the algebraic constraint in the DAE system. 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 described by the algebraic constraints. In the algorithm, the equations are never explicitly projected, but the projection is only applied as needed. Numerical experience indicates that the algorithmic choices for the control of inexactness and line-search can help avoid subproblems with matrices that are only marginally stable. 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_2020a,
    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 = {Efficient Solution of  Large-Scale Algebraic {R}iccati Equations
             Associated with Index-2 {DAE}s via the Inexact Low-Rank {N}ewton-{ADI} 
             Method},
   JOURNAL = {Appl. Numer. Math.},
  FJOURNAL = {Applied Numerical Mathematics. An IMACS Journal},
    VOLUME = {152},
      YEAR = {2020},
     PAGES = {338-354},
       DOI = {10.1016/j.apnum.2019.11.016}
       URL = {https://doi.org/10.1016/j.apnum.2019.11.016}
}

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.

@article {MHeinkenschloss_BKramer_TTakhtaganov_KWillcox_2018a,
    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},
   JOURNAL = {SIAM/ASA J. Uncertainty Quantification},
  FJOURNAL = {SIAM/ASA Journal on Uncertainty Quantification},
    VOLUME = {6},
      YEAR = {2018},
    NUMBER = {4},
     PAGES = {1395-1423},
       DOI = {10.1137/17M1160069},
       URL = {https://doi.org/10.1137/17M1160069}
}

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 show 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 = {accepted for publication}
}
	  

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 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}
}