A Direct Optimization Algorithm for Problems with Differential-Algebraic Constraints: Application to Heat and Mass Transfer

: In this article, an optimization task with nonlinear differential-algebraic equations (DAEs) is considered. As a main result, a new solution procedure is designed. The computational procedure represents the sequential optimization approach. The proposed algorithm is based on a multiple shooting parametrization method. Two main aspects of a generalized parametrization approach are analyzed in detail: a control function and DAE model parametrization. A comparison between the original and modiﬁed DAEs is made. The new algorithm is applied to solve an optimization task in heat and mass transfer engineering.


Introduction
Currently, the dynamic development of engineering systems of parametric models can be observed. The parametric models, built from variables with a known interpretation connected with the equations of the laws of physics, are able to more and more accurately reflect the studied dependencies. The development of this type of model is limited only by the scope of natural science and computing capabilities. The parallel development of the mathematical models and numerical simulation algorithms can be easily observed in many fields of engineering, in particular in chemical engineering [1] and mechanical engineering [2]. Advanced computational techniques for the process design and optimization can have a positive impact on their quality and profitability.
Literature studies indicate the development of the parametric mathematical models of heat and mass transfer phenomena, which can be applied to HVAC (heat, ventilation, and air-conditioning) systems. The design research is oriented toward the systems operating under specific environmental conditions, as well as the simulation analysis of the system performance at various operating points in order to increase the efficiency of heat and mass exchange units. Therefore, the latest achievements in the use of the mass and heat transfer models are presented.
In their work, Pandelidis et al. [3] numerically analyzed the M-cycle cooling tower operation. In this research, critical heat and mass transfer process conditions were considered. Moreover, some critical factors with the highest impact on the system performance were established. Finally, a higher than 100% wet-bulb effectiveness for water cooling was achieved.
Yang and co-workers [4] designed some theoretical models of buoyancy-driven flows by area heat sources. According to the analyzed dependencies, a new law of thermal stratified flow for the transition zone was obtained. The numerical simulation approach was used to verify the correctness of the approximate theoretical solutions.
in Section 5. The results of numerical simulations are presented in Section 6. Finally, the conclusions are given in Section 7.

Optimization with Differential-Algebraic Constraints
The starting point for further considerations is the statement of a dynamic optimization task with the differential-algebraic constraints. In some practical applications, the DAE system can be presented in a semi-explicit form:ẏ (t) = f (y(t), z(t), u(t), p, t) 0 = g(y(t), z(t), u(t), p, t) where the specified functions and variables have appropriate properties and perform the assumed roles in the model; in particular, y(t) ∈ R n y is a vector of differential state trajectories, andẏ(t) = dy dt . Moreover, z(t) ∈ R n z is a vector of algebraic state trajectories. u(t) ∈ R n u denotes the vector of control functions. p ∈ R n p is a vector of parameters constant in time. t ∈ R is an independent variable with t ∈ [t 0 t f ]. The important assumptions are connected with the functions f and g and the vector of consistent initial conditions [y(t 0 ) z(t 0 )] T = [y 0 z 0 ] T . Assumption 1. The initial conditions [y(t 0 ) z(t 0 )] T = [y 0 z 0 ] T for the system (1) are consistent, iff for a given initial value of the control function u(t 0 ), the equality: is fulfilled.
The consistent initialization of the differential-algebraic systems is an important and classical task considered in the literature. The vector [y(t 0 ) z(t 0 )] T can be obtained as a result of analytical considerations [12], characterized by the real-life process initial conditions [13], computed by the deterministic optimization procedure [14] or obtained by heuristic optimization algorithms [15].
The relation between the initial conditions' consistency and the index of the considered system can be obtained.

Assumption 2.
Let us assume that the following C 2 class functions f and g are considered: f : R n y × R n z × R n u × R n p × R → R n y g : R n y × R n z × R n u × R n p × R → R n z (3) Theorem 1. Iff a matrix ∂g ∂z is nonsingular, then the system (1) is an implicit ODE of the following form: (4) and has index-1.
Proof of Theorem 1. Let us differentiate the algebraic part of Equation (1) with respect to t: which results in: and: This is true, iff a matrix ∂g ∂z is nonsingular.
The minimum number of times that all or part of the system (1) must be differentiated with respect to t in order to determineż(t) as a continuous function of z and t is the index of the DAE (1) [16].
From the computational methods' point of view, the consistency of the proposed initial conditions and the index of a system are the main properties of the differential-algebraic equations. The index-1 assumption is common for many real-life applications, although some index-2, as well as higher index formulations are the subject of research. Throughout this work, it is assumed that the considered model constraints can be described by the system of index-1 highly nonlinear DAEs in a semi-explicit form (1).
In the context of the differential-algebraic model constraints, a scalar-valued performance index in a general form can be considered: and in particular: To optimize the performance index (8) subject to the DAE system (1), both indirect and direct methods can be applied. The indirect methods are based on the solution of the necessary optimality conditions. Thus, for the indirect methods, it is necessary to explicitly derive a system of the Euler-Lagrange equations consisting of adjoint equations, control equations, and transversality conditions, separately for each considered task [17]. Although the indirect methods have some known disadvantages, the recent research results indicate efficient indirect schemes [18], as well as combinations of direct-indirect approaches [19].
The optimization of systems described by the nonlinear differential-algebraic constraints requires a reliable approach to obtain the state trajectories for the given consistent initial conditions. Therefore, the further considerations related to the computational optimization with the nonlinear DAE model are based on the multiple shooting approach. The application of the multiple shooting step needs to divide the independent variable domain into an assumed number N of subintervals: according to the following rule: (11) and can be used to parametrize the dynamic optimization task. An appropriate optimization problem parametrization enables us to apply efficient nonlinear optimization procedures and obtain a suboptimal solution even for highly nonlinear and possibly unstable systems. There are three main dynamic optimization task approaches 1.
Control vector parametrization (CVP): The parametrization is related only to the control function and results in a medium-scale nonlinear optimization task [20].

2.
Direct shooting method: The parametrized control function, as well as the initial conditions of the state trajectories are treated as additional decision variables. The differential-algebraic equations are solved for each subinterval independently. Therefore, an efficient DAE solver is required. Additionally, the equality constraints providing the state trajectory continuity need to be taken into account [21].

3.
Direct transcription formulation: The values of the state trajectories are parametrized and treated as additional decision variables. To parametrize the state trajectories, the Runge-Kutta, Hermite-Simpson, or any other numerical scheme can be used. The direct transcription formulation results in a large-scale and sparse nonlinear optimization task. Moreover, only the final result can be considered as an applicable solution [17].
Recently, a new family of methods was constituted by a combination of the direct shooting method with a DAE model approximation. The multiple shooting approach enables us to consider other systems in each subinterval. Moreover, the appropriate approximation makes the highly nonlinear DAE models possible to solve. Then, the new model properties can be analyzed from a computational optimization point of view.

A Generalized Parametrization Approach
In the presented work, the multiple shooting approach is used to parametrize both the control function, as well as the differential-algebraic model of the system. Usually, in the direct optimization algorithms, a partial state trajectory parametrization is often used. Therefore, to unify parametrization issues, the general method for the model, state, and control function parametrization is defined. The presented discussion is an extension for the understanding of the direct shooting approach.

Definition 2. Let us denote by
, such that f : R n x → R n f and P ∈ R n f ×m is a matrix of the function f parameters, where: for x ∈ [x 0 x f ] and according to the chosen norm. In particular, the parametrization of a scalar function f (x) such that f : R n x → R is: Then, according to Definition 2, the control function u(t) ∈ R n u can take a new form: and briefly: The parameter matrix P u should be found with the condition: for t ∈ [t 0 t f ] and according to the chosen norm.
Moreover, the parametrization procedure P can be applied to the DAE constraints (1). In other words, it can be used with regard to the model relations, not the state trajectories. The parametrization procedure P can result in a new form of the model dependencies: . . .
. . . where: Finally, in the considered DAE model (17), the parametrized control function u(P u ; t) can be introduced: Let us consider the special case, when P m = P n y +n z +n u . Then, the DAE model (1) takes a special parametrized form: . . . where: · · · a n y +n z +n u 1 a 1 2 a 2 2 · · · a n y +n z +n u 2 n y a 2 n y · · · a n y +n z +n u n y a 1 n y +1 a 2 n y +1 · · · a n y +n z +n u n y +1 a 1 n y +2 a 2 n y +2 · · · a n y +n z +n u n y +2 n y +n z a 2 n y +n z · · · a n y +n z +n u n y +n z with the condition: As was indicated, the parametrization P n y +n z +n u is applied to the model, not to the state trajectory. The appropriate size of the matrix P DAE can be used to obtain a new, linear model of the DAE constraints: and similar to the state space equations. The matrix P DAE (21) manipulations can result in other forms of the considered dynamical system. Finally, the control vector parametrization P m u can be formally applied in the parametrized model (23) where A, B, and P u represent matrices of a known size, but with unknown values of their components.

The New Computational Procedure
The presented, general approach to the differential-algebraic constraint parametrization is used to design a new computational procedure. The described algorithm is not problem-specialized. Therefore, the appropriate parametrization approach predestines it to be used in various areas of technology.
The formulation of the new solution procedure is based on three aspects of the presented considerations: • the multiple shooting approach, • the DAE constraint parametrization, • the similarity between the original DAE and the parametrized system.
There are compatible mile-stones, representing appropriate steps in the new optimization task formulation. All presented aspects influence the final form of the modified optimization task, which should be characterized by a high stability and make a nonlinear dynamics system simulation possible. Therefore, it will be clearly indicated how the multiple shooting-based parametrization can be used to obtain the modified optimization problem with features similar to the original one.

The Multiple Shooting Approach
The multiple shooting approach MSM N enables us to divide the independent variable domain into an assumed number of N subintervals. Then, the relation between obtained subintervals takes the form of Equations (10) and (11). Moreover, in a sense, the differential-system of Equation (1) is divided into subintervals also. In the other words, the DAE systems can be considered independently in subintervals. This situation is presented with details in Equation (25): The new system obtained by the MSM N approach consists of the following components: • i = 1,. . . , N is the number of the current subintervals. • The current i-th subinterval is considered in a time interval (TI) The initial conditions (ICs) of the differential state trajectories are treated as decision variables and take a structured concatenated form: • EQsdenotes the system of the DAE constraints connected to the i-th subinterval. The considered system consists of the vector-valued functions f i and g i such that: and y i (t i ) ∈ R n y i is the vector of the differential state trajectories in the i-th subinterval. z i (t i ) ∈ R n z i is the vector of the algebraic state trajectories in the i-th subinterval. u i (t i ) ∈ R n u i denotes the vector of the control functions in the i-th subinterval. p ∈ R n p is the vector of the parameters constant in time. t i ∈ R is the independent variable, and It is assumed that: and this is true, e.g., for DAE models in multi-stage systems [22]. • To solve the DAE system on the appropriate subintervals (25), the consistent initial conditions (CICs) are needed. In particular, the consistent initial conditions can be obtained with the following equations:  • The continuity conditions: The multiple shooting-based algorithm should result in continuous trajectories of the state and control variables, unless explicitly stated otherwise. This means that the value of the state trajectory at the end of the current subinterval is equal to the state trajectory at the beginning of the next subinterval. The indicated relation takes the form of an equality constraint vector: 

The DAE Constraints' Parametrization
The presented MSM N -based step of the new optimization task formulation does not indicate the exact way of the state and control trajectories' parametrization. These aspects will be considered in this section.
The parametrization used should take into account the properties of the considered optimization task, in particular: • the physical (technological) properties of the system under consideration (e.g., oscillations, discontinuities, a nonlinear behavior like hysteresis, regions of instability), • the properties of the DAE solver used (the ability to determine the consistent initial conditions, the model solution accuracy in the instability region), • the performance of the optimization algorithm (the available operating memory, the possibility of parallel computations, dedicated numerical algebra procedures for sparse matrices and large-scale tasks, the globalization strategy, the rate of convergence).
The solutions presented in the literature represent two possible ways to parametrize the differential-algebraic constraints system (1).

1.
Parametrize both the state and control trajectories, but without parametrization of the relations represented by the f and g functions. This approach results in a large-scale algebraic system of equations. Let us denote: where m y , m z , m u ∈ N + and P y , P z , P u are the matrices of the approximating polynomial coefficients. Then, if:˙ y(P y ; t i ) = d y(P y ; t) dt , with the relation t 0 < · · · < t m , where m = max{m y , m z , m u }, the following algebraic system of equations G(P y , P z , P u ) is obtained: . . . y(P y ; t m ) -f ( y(P y ; t m ), z(P z ; t m ), u(P u ; t m ), p, t m ) g( y(P y ; t 0 ), z(P z ; t 0 ), u(P u , t 0 ), p, t 0 ) . . . g( y(P y ; t m ), z(P z ; t m ), u(P u ; t m ), p, t m ) The parametrization of the constraints represented by the functions f and g: Opposite the previous point, the state and control trajectories can be calculated by a numerical procedure implemented in a DAE solver. Recently, a special case of DAE model parametrization was used to introduce the variability constraints into the optimization task [22].
The choice of a particular parametrization method is dependent on the task under consideration, as well as the efficiency of the solution method used. The situation when the DAE constraints cannot be solved with the available implemented numerical procedure seems to be particularly interesting. Further considerations will be made to resolve this issue.
Let us consider the nonlinear system of differential-algebraic constraints (1). In computer simulations of many real-life problems, the direct application of the numerical DAE solvers cannot result in an accurate and stable solution. These cases were reported, e.g., in the monograph [22]. Therefore, it is suggested to apply a simplified DAE model, which can reflect the dynamical features appropriately and does not exhibit numerical difficulties or instabilities. The commonly used methodologies in the context of a highly nonlinear DAE constraint solution are based on a polynomial representation of the state trajectories y(t) and z(t) ( [23]) or a full state parametrization (direct transcription method [17]). There is a common way to transform a difficult problem of solving the DAE system into a system of nonlinear algebraic relations. Then, the iterative methods can be used to improve the current solution step-by-step [24].
Opposite the presented procedures, the direct shooting approach (a combination of the multiple shooting method with the DAE solver) can be applied. It is assumed that the solution of the DAE system obtained on a shorter subinterval can be more accurate than a solution calculated for a long domain. However, even in this reduced situation, the DAE solver can fail and result in no solution, and the optimization procedure can be interrupted. Finally, the number of shooting points should be increased, and the whole solution procedure needs to be started again. Therefore, a combination of the multiple shooting method with the modified (parametrized) DAE model solved by an outer numerical procedure is considered.
In particular, the modification of the DAE system is not only understood as its linearization: where t λ ∈ [t 0 t f ]. The linearization of the model constraints should be performed numerically or analytically at the beginning of the computations. Moreover, the obtained linearized model is dependent on the chosen simulation point (y(t λ ), z(t λ ), u(t λ )). One of the main ideas considered in this research is to estimate the unknown model parameters. This approach enables us to apply all forms of mathematical models, which can be effectively treated by the implemented outer DAE solver. The proposed parametrization approach results in the following theorem.
Theorem 2. The parametrization of the DAE system in a semi-explicit form (1) with the P m f and P m g parametrizations for the differential and algebraic model constraints, respectively, introduces: decision variables, where m f , m g ∈ N + .
Proof of Theorem 2. Let us note that: = f (P f ; y(t), z(t), u(t), p, t) g(P g ; y(t), z(t), u(t), p, t) , where m f , m g ∈ N + and P f ∈ R n y ×m f , P g ∈ R n z ×m g are the matrices of unknown parameters.

The Similarity between the Original and Parametrized DAE Systems
An important issue is related to the method of measuring the similarity between the original DAE system and the modified model with the estimated coefficient matrices P f and P g . As was mentioned in the previous section, the difference between the systems can have an integral form adapted to the multiple shooting method, according to the following definition. Definition 3. The difference between the original and modified DAE systems in a semi-explicit form (1) can be expressed as: where · denotes a chosen norm, i is the index of a subinterval, and j denotes the index of the function vector component.
Further considerations require the following assumption.
Assumption 3. Functions f i j (·, t i ), f i j (·, t i ), g i k (·, t i ) and g i k (·, t i ) for i = 1, . . . , N, j = 1, . . . n y , k = 1, . . . , n z are continuous, real-valued, and defined on the appropriate subintervals t i ∈ [t i 0 t i f ] for all considered consistent initial conditions and control functions u i (t i ) and u i (·, t i ).
Based on Assumption 3, let us consider the two following DAE constraint systems: and:˙y (t) = f (P f ;y(t),z(t), u(P u ; t), p, t) 0 = g(P g ;y(t),z(t), u(P u ; t), p, t) To estimate the difference between the differential state trajectories: the implications of the mean value theorem can be applied.
Theorem 3. Ifẏ(t) =ẏ(t) for all t in an interval t ∈ (t 0 t f ) of the domain of these functions, then y(t) −y(t) is constant or y(t) =y(t) + c, where c ∈ R is constant on (t 0 t f ).
As a result of Theorem 3, to supervise the fulfillment of the condition: the vector of equality constraints is considered as a part of the new optimization task: The values of z i (t i 0 ), i = 1, . . . , N are a consequence of the consistent initial conditions.

The New Solution Procedure
The main steps of the new solution algorithm are described briefly as "PRE Procedure −−−−−→ POST", where "PRE" denotes an input of the "Procedure". Then, "POST" is the result of the "Procedure" for the input "PRE".
The parametrized DAE-multiple shooting method algorithm: 1. Define the index-1 DAE constraints in the form (1). 2. For the system (1), find the values n y , n z , n u ∈ N + . 3. Define the number N ∈ N + , which denotes the number of considered subintervals: 4. Implement the continuity constraints in the form (34). 5. Perform the control vector parametrization. Define m u ∈ N + , then: where P u ∈ R (n y +n z )×m u . 6. The DAE model parametrization: Define m f , m g ∈ N + , then: where P f ∈ R n y ×m f and P g ∈ R n z ×m g . 7. The DAE models' similarity: Define the vector of the equality constraints according to Equation (46).
and, in particular: where t i ∈ [t i 0 t i f ] and i = 1, . . . , N. 9. The nonlinear optimization task: The objective function (50) is optimized subject to the modified DAE constraints (48) and (49). Then, the vector of the consisted initial conditions is needed to perform the numerical computations correctly. Additionally, the continuity and similarity constraints are considered. The new nonlinear optimization task is formulated by the indicated components. The numerical continuous optimization procedure, which co-operates with the outer DAE solver, can be used to design an efficient sequential optimization algorithm. 10. The trajectories y(t), z(t), and u(t), where: as well as the appropriate matrices P f i , P g i , P u i and the vector of the initial conditions x y i for i = 1, . . . , N are the solutions obtained by the new optimization algorithm. Moreover, it is suggested to save the final value of the similarity constraint vector. This information can be used to improve the similarity between the original and modified DAE constraints by a modification of the subinterval number N.

An Example of the Numerical Calculations
The designed algorithm is applied to optimize the heat and mass transfer phenomena described by a system of nonlinear differential-algebraic equations. In general, the technological effectiveness of the presented system in comparison with other available configurations was indicated in [25]. The modified counter-flow exchanger design task [26] with the DAE constraints in a semi-explicit form is constructed as follows: min subject to:ẏ 1 (l) = −B · (z 1 (l) − y 1 (l))/l ḟ y 2 (l) = C · (z 2 (l) − y 2 (l))/l ḟ y 3 (l) = C · (z 3 (l) − y 3 (l))/l f 0 = E 1 ·ẏ 1 (l) −ẏ 2 (l) − E 2ẏ3 (l) 0 = z 2 (l) − (z 1 (l) − D · (y 1 (l) − z 1 (l))) 0 = z 3 (l) − 0.622 · z 4 (l) P b −z 4 (l) ) 0 = z 4 (l) − 6.107 · e 0.0726·z 2 (l)−2.912·10 −4 ·(z 2 (l)) 2 +8.33·10 −7 ·(z 2 (l)) 3 where y 1 (l)-temperature in the first channel ( • C), y 2 (l)-temperature in the second channel ( • C), y 3 (l)-humidity ratio in the second channel (kg/kg), z 1 (l)-temperature of the plate surface in the first channel ( • C) z 2 (l)-temperature of the plate surface in the second channel, z 3 (l)-humidity ratio on the plate surface in the second channel (kg/kg), and z 4 (l)-the static pressure (h Pa). Moreover, there is a vector of global parameters: This model consisting of the differential-algebraic constraints has nonlinear algebraic relations. Therefore, analytical approaches are difficult to apply directly. Moreover, the numerical difficulty depends on the assumed parameters in the vector p. In particular, the values of B and C show how dynamic the mass and heat transfer is in the system. Then, according to the vector p, the solution procedure can result in an accurate solution for a slow dynamics or it can result in a computing error, if a considered transfer is rapid or sharp. The performed calculations indicate that the single shooting approach can be used if B and C are less than 10. The direct shooting strategies can be useful for B and C less than 20. For more dynamic processes, the co-operation of the multiple shooting with model parametrization seems to be appropriate.
To solve the presented optimization task, the multiple shooting approach with DAE constraint parametrization is implemented. The applied parametrization fulfills the conditions: and: The algorithm is implemented with two numerical procedures from the MATLAB environment [27]: • fmincon: The numerical optimization procedure is used to solve the parametrized continuous optimization task subject to the equality constraints. The numerical values of the objective and constraint functions are obtained with the outer DAE solver. • ode15s: The outer solver is used to supply the state trajectories of the modified DAE relations for the given consistent initial conditions and values of parameters. This solver gives a support to finding consistent initial conditions.
The new algorithm with 10 equidistant subintervals results in a stable solution. The added constraints force the continuity of the state trajectories, which are presented in Figures 1-3. The solution obtained for each subinterval is marked with a different color. The initial value of y 1 (0) = 17.7556 • C is obtained. The value of the objective function is sequentially minimized, and finally, (y 1 (l f ) − 30) 2 ≤ 10 −6 .
As was indicated previously, for the considered task with the parameter vector p (57), the classical single shooting-based approach cannot obtain any useful solution. Therefore, the efficient modifications of the multiple shooting method enable us to have a detailed insight into the course of the phenomena in a wide range of model parameters. The simulations can be performed even for such conditions, which seems to be difficult from the numerical calculations' perspective.

Conclusions
In the presented research, a new algorithm for the optimization task subject to the nonlinear differential-algebraic constraints is designed. The introduced methodology is an extension of the sequential approach for optimization, where the outer DAE solver can be used and the obtained state and control trajectories can be modified by the numerical optimization procedure step-by-step.
As a main result of this work, the new nonlinear optimization task formulation is derived. This formulation can be used and further analyzed in the context of the modern numerical optimization procedures with equality constraints. In other words, all parts of the nonlinear optimization task are clearly indicated and can be further modified and developed.
Much attention has been paid to the method of limiting the system dynamics in order to enable its numerical simulation. Therefore, the proposed method is based on a parametrization of the considered DAE constraints. The generalized parametrization approach is discussed in detail. Because of the system parametrization, during the optimization procedure, the similarity between the original and modified (parametrized) DAE constraints is measured. Moreover, the similarity constraints are considered as a part of the nonlinear optimization task.
The physical model parameters of heat and mass transfer systems have clear and known interpretations. In this work, a solution is sought for a such problem, which so far could not be solved for certain parameter values. Now, the values of the vector p can be modified according to a wide range of industrial applications. The simulation studies will enable the design of new mass and heat exchangers. In the presented case study, the solution is obtained with the co-operation of the numerical procedures fmincon and ode15s in the MATLAB environment. However, in fact, the choice of numerical methods is not imposed, but may be based on available procedures.
Some of the highlighted issues will be detailed in the future. Further work will be concentrated on issues related to special forms of the DAE model parametrization, as well as their properties and relations with the original dynamics. The extended similarity testing procedures will be considered also. The presented algorithm can find successful applications in other fields of engineering, like the design of rapid thermal processes in aviation, space, and nuclear energy.