Simulation of differential-algebraic equation systems with optimization criteria embedded in Modelica

https://doi.org/10.1016/j.compchemeng.2020.106920Get rights and content

Abstract

Differential-algebraic equations with embedded optimization criteria (DAEO) are a class of mathematical models for underdetermined differential-algebraic equation (DAE) systems with less algebraic equations than algebraic variables. The algebraic variables may be calculated as the solution of an embedded (non)linear program, yielding a DAEO system. An example for DAEOs is the dynamic flux balance analysis (DFBA) approach, where the formulation of metabolic reaction networks leads to an underdetermined equation system for the intracellular fluxes that are assumed to behave optimally with respect to some cell-specific optimization criterion.

We present a toolbox that allows formulation of DAEOs in the object-oriented Modelica modeling language. The solution method is based on substituting the embedded optimization problem with its first-order Karush-Kuhn-Tucker conditions to obtain a nonsmooth DAE system that can be simulated by a root-finding DAE solver. One nonlinear example and two examples based on DFBA demonstrate the performance of the toolbox.

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.,y˙d(t)=f(yd(t),ya(t),t,p),withyd(0)=yd,0(p)ya(t)argminyaRnah(yd(t),ya,t,p)s.t.0=gk(yd(t),ya,t,p),k=1,,ne0gk(yd(t),ya,t,p),k=ne+1,,ng

Here, yd(t)Rnd and ya(t)Rna are the differential and algebraic variables, respectively. The differential equations are given by sufficiently smooth functions f:Rnd×Rna×[0,tf]×RnpRnd. The vector pRnp 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., h=cTya 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)

  • K.E. Brenan et al.

    Numerical solution of initial-value problems in differential-algebraic equations

    (1996)
  • S.L. Campbell et al.

    Solvability of general differential algebraic equations

    SIAM Journal on Scientific Computing

    (1995)
  • A. Caspari et al.

    Dynamic optimization with complementarity constraints: Regularization for direct shooting

    Optimization Online

    (2019)
  • F.E. Cellier et al.

    Continuous System Simulation

    (2006)
  • N.P. Faisca et al.

    Multiparametric Linear and Quadratic Programming

    (2011)
  • H. Ferreau et al.

    An online active set strategy to overcome the limitations of explicit mpc

    International Journal of Robust and Nonlinear Control

    (2008)
  • A.V. Fiacco

    Introduction to sensitivity and stability analysis in nonlinear programming

    (1983)
  • P. Fritzson et al.

    Modelica - a general object-oriented language for continuous and discrete-event system modeling and simulation

    Proceedings 35th Annual Simulation Symposium

    (2002)
  • J.A. Gomez et al.

    DFBAlab: A fast and reliable MATLAB code for dynamic flux balance analysis

    BMC bioinformatics

    (2014)
  • V. Gopal et al.

    Smoothing methods for complementarity problems in process engineering

    AIChE Journal

    (1999)
  • J. Guddat et al.

    Parametric Optimization: Singularities, Pathfollowing and Jumps

    (1990)
  • T.J. Hanly et al.

    Dynamic flux balance modeling of microbial co-cultures for efficient batch fermentation of glucose and xylose mixtures

    Biotechnology and bioengineering

    (2011)
  • R. Hannemann-Tamás et al.

    Adjoint sensitivity analysis for nonsmooth differential-algebraic equation systems

    SIAM Journal on Scientific Computing

    (2015)
  • S.M. Harwood et al.

    Efficient solution of ordinary differential equations with a parametric lexicographic linear program embedded

    Numerische Mathematik

    (2016)
  • J.L. Hjersted et al.

    Optimization of fed-batch Saccharomyces cerevisiae fermentation using dynamic flux balance models

    Biotechnology progress

    (2006)
  • Cited by (6)

    • Nonlinear programming reformulation of dynamic flux balance analysis models

      2023, Computers and Chemical Engineering
      Citation 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 Engineering
    • Direct single shooting for dynamic optimization of differential-algebraic equation systems with optimization criteria embedded

      2022, Computers and Chemical Engineering
      Citation 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.

    View full text