Simulation of differential-algebraic equation systems with optimization criteria embedded in Modelica
Introduction
In chemical engineering, dynamic models are important for an improved understanding of transient process behavior, in particular in response to disturbances. They also help to design startup and shutdown procedures and form the basis for model-based operation and control. In some cases, however, the system of differential-algebraic equations (DAE) may be underdetermined due to for instance lack of mechanistic knowledge, i.e., there may be more algebraic states than algebraic equations.
The simulation of such underdetermined DAEs is a non-trivial problem. One may interpret these systems as differential inclusions where the right-hand side of the differential equation is a set-valued function (Smirnov, 2002). Under some assumptions (e.g., linear differential and algebraic equations), it might be possible to obtain a closed-form description of the set of trajectories. However, the solution of general nonlinear underdetermined DAEs is much more challenging. Methods based on Monte Carlo sampling or methods that approximate the reachable set (Villanueva et al., 2014) have to be employed.
Nevertheless, from the modeling perspective it is often favorable to impose one optimization condition which selects a unique or at least locally unique solution from the space of possible solutions. This approach leads to a differential-algebraic equation system with optimization criteria embedded (DAEO). This problem class is relevant to a variety of chemical engineering applications.
For instance, in systems biology, a successful method to model the complex reactions inside a living organism is based on steady-state balance equations coupled with the stoichiometry of intracellular reactions and additional thermodynamic, kinetic or phenomenological constraints to narrow the possible metabolic flux distributions inside the cell. In most cases, the resulting (linear) equation system is underdetermined. The flux balance analysis (FBA) method selects one possible flux distribution by imposing an optimization criterion subject to these constraints Varma, Palsson, 1994, Palsson, 2008. The dynamic flux balance analysis (DFBA) method couples the FBA approach to the dynamics of the reaction medium Mahadevan, Edwards, Doyle, 2002, Palsson, 2008. DFBA became one of the dominating methodologies used to describe the transient behavior of large biochemical networks inside living cells, for example, to model the growth of Escherichia coli (E. coli) on d-glucose and acetate (Mahadevan et al., 2002) or to investigate growth and ethanol productivity of a co-culture of E. coli and Saccharomyces cerevisiae on d-glucose/d-xylose mixtures (Hanly and Henson, 2011). The DFBA approach leads to a DAEO where the equality and inequality constraints of the embedded optimization problem are linear. The objective function may be linear (e.g., maximization of cell growth Varma, Palsson, 1994, Hanly, Henson, 2011 or nonlinear (e.g., maximization of biomass yield per flux unit (Zhao et al., 2017).
Dynamic models of separation processes often assume thermodynamic equilibrium between coexisting phases (e.g., vapor-liquid equilibrium) that is found at the (global) minimum of Gibbs free energy Baker, Pierce, Luks, 1982, Michelsen, 1982, Michelsen, Mollerup, 2007. Here, the phase or chemical equilibrium assumption is made on purpose due to much faster mass exchange dynamics compared to the overall dynamics of the system. In contrast, the embedded optimization problem in the DFBA approach is motivated by a lack of knowledge regarding the intracellular flux distribution. Nevertheless, dynamic models based on phase equilibrium may be formulated as DAEOs where both objective function and constraints (e.g., material balance) are nonlinear. The applications range from modeling of organic aerosol particles (Landry et al., 2009) to dynamic simulation of separation processes based on different types of phase equilibria Gopal, Biegler, 1999, Biegler, 2010, Sahlodin, Watson, Barton, 2016.
There are different approaches to solve DAEOs. The direct approach (DA) is commonly used if an embedded LP is considered Hanly, Henson, 2011, Hjersted, Henson, 2006. It calculates the algebraic variables by directly solving the optimization problem in each time step of the integration. The direct approach is used in some MATLAB implementations, for example, DyMMM (Zhuang et al., 2011) and as part of the COBRA toolbox (Schellenberger et al., 2011). While this approach is computationally fast enough for a certain class of DAEOs (in particular embedded LPs), it bears the risk of wrong simulation results caused by the interaction of optimization and integration algorithm. For example, an attempted integration step may lead to infeasibility of the embedded optimization problem (Harwood et al., 2016). A remedy is to use first-order optimality conditions (KKT conditions) to reformulate the DAEO into a nonsmooth DAE. This idea has been used to solve DAEO applications with NLP embedded (Landry et al., 2009) and resulted in Fortran (Höffner et al., 2013) and MATLAB (Gomez et al., 2014) implementations for solving models with LPs embedded.
We note that for efficient simulation of DAEOs taylored solution methods are required. These solution methods should exploit the fact that during numerical integration of the models, a sequence of related linear, quadratic or nonlinear programs has to be solved. Actually, we have parametric optimization problems, with the time as varying parameter. A comprehensive overview about the properties of parametric nonlinear programs is given by Guddat et al. (1990). If the embedded optimization problem is linear or quadratic, efficient solution methods could borrow ideas from discrete time linear model predictive control, e.g., from multi-parametric programming (Faisca et al., 2011) or online active set strategies for QP solvers (Ferreau et al., 2008). The latter will find a resemblance in our subsequently presented active-set strategy for solving DAEOs with linear or quadratic programs embedded. However, it is important to stress the difference between DAEOs and model-predictive control. In the former, an embedded optimization problem has to be solved, in the latter we have an outer optimization problem. The similarity arises because both are in fact parametric optimization problems.
In comparison to existing literature, this work does not focus on a particular application but aims to treat general DAEOs. Thereby, we build on the solution technique proposed by Landry et al. (2009) using the first-order KKT conditions to reformulate general DAEOs into nonsmooth DAEs and provide implementations of the direct approach and the KKT embedding approach. We build on Modelica as an object-oriented modeling language that provides a structured and efficient approach to model systems with differential, algebraic and discrete equations. For solving the embedded optimization problem, we implemented an external toolbox that is interfaced to the Modelica-based simulation environment. Our point of view is that Modelica is better suited to model complex systems than MATLAB or Fortran. The aim of Modelica is to unify and generalize previous object-oriented modeling languages and become a de facto standard (Fritzson and Bunus, 2002). Cellier and Kofman compare Modelica and MATLAB for the simulation of an electrical circuit and conclude that in MATLAB much more manual preprocessing is required (Cellier and Kofman, 2006). From our personal experience, setting up a simulation with Fortran libraries is even more tedious.
By providing a DAEO toolbox for Modelica, we get an easy-to-use interface for model extension and embedding of DAEO formulations into large and complex models (e.g., DFBA-based bioreactor coupled with plant-wide biorefinery model (Ploch et al., 2019)).
The structure of this work is as follows. In the next section, we formally introduce the DAEO formulation, describe the applied solution methods and discuss the specialized solution method for DFBA models and the implementation within the Modelica language. Subsequently, we demonstrate the applicability of our method on one nonlinear example and two linearly constrained examples based on DFBA. Finally, we close with concluding remarks. The DAEO software is available on the public git repository with the permalink http://permalink.avt.rwth-aachen.de/?id=252659.
Section snippets
DAEO formulation
In DAEOs, the algebraic variables are determined by an embedded nonlinear program (NLP), i.e.,
Here, and are the differential and algebraic variables, respectively. The differential equations are given by sufficiently smooth functions . The vector comprises all model parameters, tf is the final time and t ∈ [0, t
Specialized solution method for DFBA models with linear objective
In the context of DFBA models, the constraints of the embedded optimization problem (1c) and (1d) are linear and only a small part of the formulation is time-dependent. In this section, we specialize the KKT embedding solution technique for this modeling class to enable fast and reliable simulation. In particular, we propose a regularization to avoid non-unique solutions and exploit the structure of the model formulation. In case of a linear objective function (i.e., with objective
Illustrative examples
The following examples demonstrate the modeling of DAEOs and their solution using the nonsmooth system approach described in Section 2.3 and its advantages in terms of required CPU time and the possibility of provide analytic derivatives (using the implicit function theorem) compared to the direct approach (cf. Section 2.2). The first example comprises the numerical results for the small nonlinear DAEO (2). The following two examples are based on the DFBA approach and discuss the different DAEO
Conclusion and outlook
This work presents and compares different solution methods for differential-algebraic equation systems with embedded optimization criterion. The implemented DAEO toolbox allows model formulation in the modeling language Modelica. The direct solution method directly optimizes the embedded (non)linear program in each step and successfully solves the considered examples. The KKT embedding approach reformulates the DAEO into a nonsmooth DAE by using first-order optimality conditions. The
Acknowledgements
This work was supported by the German Research Foundation (DFG) under grant MI 1851/3-1.
References (50)
The isothermal flash problem. Part I. Stability
Fluid Phase Equilibria
(1982)- et al.
Numerical optimization
(2006) On an existence and uniqueness theory for nonlinear differential-algebraic equations
Circuits, Systems, and Signal Processing
(1991)- et al.
Quantitative prediction of cellular metabolism with constraint-based models: The COBRA toolbox v2.0
Nature protocols
(2011) - et al.
Thermodynamics and its applications
(1997) - et al.
On the stability of set-valued integration for parametric nonlinear ODEs
24th European Symposium on Computer Aided Process Engineering
(2014) - Zhao, X., Ploch, T., Mitsos, A., Noack, S., Wiechert, W., von Lieres, E., 2020. Analyze the local well-posedness of...
- et al.
Introduction to numerical continuation methods
(2003) - et al.
Gibbs energy analysis of phase equilibria
Society of Petroleum Engineers Journal
(1982) Nonlinear programming: Concepts, algorithms, and applications to chemical processes
(2010)
Numerical solution of initial-value problems in differential-algebraic equations
Solvability of general differential algebraic equations
SIAM Journal on Scientific Computing
Dynamic optimization with complementarity constraints: Regularization for direct shooting
Optimization Online
Continuous System Simulation
Multiparametric Linear and Quadratic Programming
An online active set strategy to overcome the limitations of explicit mpc
International Journal of Robust and Nonlinear Control
Introduction to sensitivity and stability analysis in nonlinear programming
Modelica - a general object-oriented language for continuous and discrete-event system modeling and simulation
Proceedings 35th Annual Simulation Symposium
DFBAlab: A fast and reliable MATLAB code for dynamic flux balance analysis
BMC bioinformatics
Smoothing methods for complementarity problems in process engineering
AIChE Journal
Parametric Optimization: Singularities, Pathfollowing and Jumps
Dynamic flux balance modeling of microbial co-cultures for efficient batch fermentation of glucose and xylose mixtures
Biotechnology and bioengineering
Adjoint sensitivity analysis for nonsmooth differential-algebraic equation systems
SIAM Journal on Scientific Computing
Efficient solution of ordinary differential equations with a parametric lexicographic linear program embedded
Numerische Mathematik
Optimization of fed-batch Saccharomyces cerevisiae fermentation using dynamic flux balance models
Biotechnology progress
Cited by (6)
Nonlinear programming reformulation of dynamic flux balance analysis models
2023, Computers and Chemical EngineeringCitation Excerpt :Another class of approaches proposes the reformulation of the dFBA using the Karush Kuhn Tucker (KKT) conditions of the FBA problem. Ploch et al. (2020) reformulated the dFBA model as a nonsmooth DAE system and applied a homotopy continuation to avoid the fluxes from jumping between optimal solutions. Scott et al. (2018) used an interior point reformulation of FBA.
Dynamic modeling of aqueous electrolyte systems in Modelica
2022, Computers and Chemical EngineeringDirect single shooting for dynamic optimization of differential-algebraic equation systems with optimization criteria embedded
2022, Computers and Chemical EngineeringCitation Excerpt :In this work, we build on the sequential approach to optimize DAEO systems. We formulate first-order optimality conditions of the embedded optimization problem to obtain a nonsmooth DAE system (Ploch et al., 2020). As previously discussed, we cannot guarantee optimality of the lower level problem.
A direct optimization algorithm for problems with differential-algebraic constraints: Application to heat and mass transfer
2020, Applied Sciences (Switzerland)