Simultaneous simulation and optimization of multiple dividing wall columns

In this work we present a new approach that we use to simulate and optimize multiple dividing wall columns at the same time. Instead of considering all model equations as constraints and all process variables as optimization variables in a large and highly nonlinear optimization problem we only incorporate a subset of the model equations as constraints and a subset of the process variables as optimization variables. The remaining process variables are calculated from this subset by a robust and fast calculation procedure. This calculation procedure also ensures that the remaining model equations are satisﬁed. A comparison with the commercial process simulator Aspen Plus R (cid:13) shows that with the new approach multiple dividing wall columns can be optimized more stable and better solutions are found. Moreover the time needed to ﬁnd an optimal design decreases signiﬁcantly.


Introduction
In process design, the problem of finding a steady state can be formulated as a system of nonlinear equations: where x var = (x j ) j∈Ivar are the process variables and g i (x var ), i ∈ I MESH are the model equations. In the case of distillation processes they are given by the socalled MESH equations (material-balance, equilibrium-balance, summation-law and heat-balance). Typically there are less equations than variables. The user has to specify a value for each remaining degree of freedom. The simulation software then tries to determine values for all process variables by solving the underlying non-linear system of equations given in Equation (1). If not just any design is to be found, the specifications are modified until a good design is found. This can for example be achieved using an optimization routine. Many classic operating units such as distillation columns or reactors can be modeled this way. They can be combined in a flow sheet and simulated safely. Introducing an objective function, the resulting flow sheets can be optimized. However there are operating units which still cause difficulties and can not be simulated safely. In order to nevertheless find a design for such units, a lot of effort is spent. For example one can try to find a thermodynamic consistent flow sheet. Good starting values are needed to find a convergent simulation. If an optimization has to be performed, the convergent area has to be explored carefully and bounds have to be set in such a way that the optimization doesn't end up with non convergent simulations. An example for these more demanding unit operations in terms of simulation are dividing wall columns Generally dividing wall columns are intensified distillation processes. Conventionally, multi-component mixtures are separated in sequences of simple distillation columns where the number of columns is one below the number of components in the original mixture. Figure 1a shows an example for the separation of a quaternary mixture using a conventional split sequence. In the figure the placeholders A, B, C and D denote the components in the feed mixture sorted according to their boiling temperature where A is the light boiler. Applying internal partition walls enables the separation of the same mixtures in only one column shell Kiss and Bildea 2011;Halvorsen et al. 2013;Jiang et al. 2018. This way pure products can also be obtained from side flows. In Figure 1b an example with three dividing walls to separate quaternary mixtures is shown. A distillation column with one dividing wall (to split ternary mixtures) is called simple dividing wall column. If more than one dividing wall is present, the column is called multiple dividing wall column (mDWC). In both cases the internal vapor and liquid flows have to be split at the partition walls. Like this, remixing of intermediate boiling components, which usually occurs in column sequences, is avoided resulting in energy savings. For the simple dividing wall the savings are around 30 % M. A. Schultz et al. 2002 and for the mDWC up to 50 % Dejanović et al. 2011;Dejanović et al. 2014.
A quarternary mixture can also be separated in a simplified version introduced in Figure 1c. It only requires two instead of three dividing walls, thus less vapor and liquid splits at these dividing walls have to be known and two sections are saved. Interestingly, this kind of column can be operated without additional energy compared to the version with three dividing walls in Figure 1b for a large set of systems Dejanović et al. 2014.
The above stated energy saving potential can only be totally exploited by an optimally designed column. For this purpose, different possibilities have been developed in the literature. The simplest way is as described above. First specifications are added to Equation (1) such that all degrees of freedom are fixed. A process simulator is used to solve the resulting system of equations. In an outer loop the specifications are modified until a satisfactory design is found. This is for example applied by Kurnatowski et al. 2017;Waltermann, Sibbing, and Skiborowski 2019. However, due to this outer loop the procedure requires significant calculation times.
An alternative to repeatedly solving the large nonlinear systems of equations with a process simulator is to embed the process simulation into an optimization problem. This can be done by incorporating all MESH equations (Equation (1)) as constraints whereas the process variables x var serve as optimization variables Biegler, Grossmann, and Westerberg 1997;Dowling and Biegler 2015. An objective (eg. minimal energy consumption) and further constraints (eg. purity demands) can be added directly to the optimization problem. The optimization solver finds a solution to the equations describing the simulation that minimizes the objective. The main challenge of this approach is that the resulting optimization problems are large and may be difficult to initialize and to solve.
To overcome these difficulties we follow in this work an approach described recently Seidel et al. 2020. The idea is to not select all but a subset of the process variables as optimization variables. Next, the structure of the MESH equations is exploited to calculate values for all process variables and presolve a part of the equations. This way one arrives at a comparably small optimization problem with favorable convergence properties.
The approach exhibits a large flexibility in the choice of the optimization variables and the choice of the MESH equations that are presolved. Indeed, both approaches described above can be modeled. Either by presolving all MESH equations and adding no equation to the optimization problem, or by using all process variables as optimization variables and adding all MESH equations to the optimization problem. But also other choices are possible. For example Kravanja and Grossmann 1996 solved complete operating units by an external function. Biegler and Hughes 1982 added only equations introduced by socalled tear streams to the optimization problem. It is also possible to calculate thermodynamic properties such as the vapor liquid equilibrium externally. This has for example be done by Brusis 2003;Poth 2003;Skiborowski, Wessel, and Marquardt 2014. Despite this flexibility, the goal is to choose a small subset of optimization variables. For this reason, not only the equations describing the thermodynamic properties but a majority of the model equations are presolved. This way only a small subset of the model equations and of the original variables need to be incorporated in the optimization problem which makes the problem small. In addition, the mapping that determines all process variables from the optimization variables must be evaluated multiple times. Thus, it must be possible to do this quickly and robustly. For this reason, most, but not all, of the equations describing a unit are solved by an external routine. The right choice depends on the considered operating unit. Hoffmann et al. 2017 presented a routine for a conventional distillation column. By Seidel et al. 2020 this was extended to the interconnection of multiple distillation columns. The goal of this paper is to develop a method for simple and multiple dividing wall columns.
The paper is organized as follows. In Section 2 we introduce the model describing a mDWC in Detail. We describe the used MESH equations and introduce an thermodynamically equivalent Petlyuk sequence. In the next section we describe the new solution approach. We first describe the general procedure and show how all process variables can be calculated from the chosen optimization variables. At the end of this section we introduce conditions which ensure that all values for the calculated concentrations remain positive. After a short description of the implementation details we present numerical results for different use cases in Section 5. First considering simple dividing wall columns, then for multiple dividing wall columns.

Equilibrium stage model for mDWC
In this section we introduce the MESH equations that describe a mDWC. For simplicity we consider a given feed and a constant pressure p within the column. However, the following can also be extended to a variable feed and pressure. We use the equilibrium stage model widely applied to model conventional distillation columns Biegler, Grossmann, and Westerberg 1997 and adjust it to the situation where dividing walls are present. The column is modeled with a total hight of N h theoretical stages. Additionally, we model a total condenser with dutẏ Q C < 0 and a total reboiler with dutyQ R > 0, which are not counted within the theoretical stages. Throughout the following, stages are numbered from bottom to top. We assume that a dividing wall starts and ends with a stage.
For a stage k in a conventional distillation column there is a liquid stream (with molar flow rate L k and molar concentrations x k ) flowing from it to the stage beneath and a vapor stream (with molar flow rate V k and molar concentration y k ) flowing from it to the stage above. Likewise, there is a liquid stream and a vapor stream flowing onto it. The equations used to describe their dependency are the MESH equations stated in the following for N C components: • Equilibrium-condition: in this work extended Raoult's law is used to model the vapor-liquid equilibrium. This means that for every i, where T k denotes the temperature in stage k, p S i the vapor pressure of the pure component i and γ i the activity coefficient for component i.
• Summation-law: • Enthalpy-balance: As we assume the equilibrium condition in each tray, the temperature can be calculated by using extended Raoult's law. Using this temperature the enthalpy of the liquid phase only depends on the concentration x and is denoted by l(x). Analogously the enthalpy of the vapor phase y is denoted by v(y).
For the feed stage of the column the material and the enthalpy-balance have to be adjusted accordingly. Three additional types of stages need to be considered at the dividing walls which are shown in Figure 2. A stage separated by a dividing wall is shown in Figure 2a. We will assume that there is no heat transfer through the dividing wall. This means that we can consider a separated stage as two independent stages k andk and formulate the MESH equations given above for each one of them. Analogously a stage divided by more than one wall is modeled by more than two independent stages. We next describe stages beneath and above a dividing wall. A stage directly above of a wall, shown in Figure 2b, has two instead of just one vapor stream entering, one from the left of the dividing wall and one from the right. The outgoing liquid stream is split into two streams, one going to the left of the dividing wall and one to the right. Both outgoing liquid streams will have the same concentrations. For convenience we model the two streams with two concentration vectors but we add the equations to ensure that their composition is equal. In the material-balance and the enthalpy-balance (Equations (2) and (5)) terms for the additional ingoing and outgoing streams have to be added. The equlibrium-condition and the summationlaw (Equations (3) and (4)) do not need to be modified as Equation (6) ensures that they also hold forx n . A stage directly beneath a dividing wall has two ingoing liquid streams and the outgoing vapor stream is split. A corresponding stage is shown in Figure 2c. Terms for these streams are added to the material-balance and enthalpy-balance. Again we model two concentrations for the outgoing vapor streams and add the equations Besides the description of the MESH equations for each theoretical stage, a mDWC can be described by a thermodynamic consistent flowsheet. Therefore sections of the mDWC are represented by fully thermally coupled simple columns without reboiler and condenser. Each of these columns has an ingoing vapor and outgoing liquid stream at the bottom of the column and an outgoing vapor and ingoing liquid stream at the top of the column. The feed stage corresponds to a split stage above or beneath a dividing wall. The columns representing sections of the mDWC are interconnected by vapor and liquid streams. The product streams are located at the top, or bottom of the representing columns. Such a flowsheet is also called Petlyuk sequence (see for example Triantafyllou and Smith 1993). An example for the separation of four components is shown in Figure 3 for the mDWC introduced in Figure 1c. The column sections are shown by dashed lines in Figure 3a. The corresponding columns in the thermodynamic consistent flowsheet are shown in Figure 3b. The column sections C i,j are numbered from right to left and from top to bottom. We skip the second index j if there is only one.

Embedding into an optimization problem
After introducing the model of a mDWC we now turn to the solution procedure. As described in the introduction we formulate the solution of the simulation as an optimization problem. Instead of considering all process variables x var = (x j ) j∈Ivar and all MESH equations g i , i ∈ I MESH in the optimization problem, we select a subset of the process variables as optimization variables x opt = (x j ) j∈Iopt with I opt ⊆ I var . Using values for these optimization variables, all process variables are calculated using a part of the MESH equations. This can be described as a map x opt → x var .
The calculated process variables already satisfy all MESH equations used in this calculation procedure, but a part of the MESH equations I constr ⊆ I MESH can be initially violated. They are incorporated in an optimization problem: where f denotes an objective and s k , k ∈ S eq and s l , l ∈ S ineq denote further constraints on the process such as purity specifications. A solution of the optimization problem satisfies all MESH equations, all specifications and minimizes the objective. We start with an initial guess for the variables x opt (Step 0) and calculate all process variables (Step 1). Using all process variables, the remaining MESH equations, the specifications and the objective function can be evaluated (Step 2).Based on their values the optimizer will suggest new optimization variables (Step 3). As the process variables depend on the optimization variables, the map given in Equation (8) has to be evaluated again (Step 1). This is repeated until all constraints are satisfied and optimality is reached. The steps of the algorithm can be summarized as follows: Step 0 : Choose initial values for x opt .
Step 2 : Evaluate missing equations and objective function.
Step 3 : Are all equations fulfilled and is optimality reached? Initial values for the optimization variables can for example be obtained using a shortcut method. A possibility for multiple dividing wall columns has been presented by Halvorsen andSkogestad 2003b andGrützner 2018. In these works a V min diagram is used to approximate the minimal amount of vapor needed to achieve the separation of the substances under consideration assuming an infinite number of stages. An improved approximation in the case of a finite number of stages has been developed by Ränger and Grützner 2021. However, in Section 5.2 we investigate the influence of the starting values on the result. The experiments show that the approach is very robust and often converges to an optimal solution even for a random choice of initial values. Nevertheless a good choice of initial values can improve the convergence rate.
For the algorithm, the evaluation of the map described in Equation (8) is crucial. It must be possible to evaluate this map reliably and fast on a large range of optimization variables x opt . Therefore, the optimization variables and the mapping to all process variables have to be chosen carefully. In the following Subsection 3.1, we will describe how to compute all process variables knowing the values for the streams that connect the column sections C i,j in the Petlyuk sequence discussed in Section 2 and shown in Figure 3 for the example of a column with two dividing walls. In Subsection 3.2 we will then introduce conditions under which this procedure is guaranteed to succeed. Thus, a suitable choice of optimization variables for multiple dividing wall columns are variables that describe the interconnecting streams of the corresponding Petlyuk sequence.
3.1 Calculation of all process variables using a stage-tostage procedure As described before the optimization routine will suggest values for all streams depicted in Figure 3 connecting the different column sections C i,j . The theoretical stages corresponding to one of these column sections is shown in more detail in Figure 4. To simplify the notation in the following, we number the stages in this column section again from 1 to n despite the stage index they already have in the complete column. The variables used as optimization variables describe the two streams at the bottom, at the top and two streams at the split or feed stage n s (shown in Figure 4 as thick lines). The type of the streams at stage n s vary between the different column sections. In Figure 4 the stage n s corresponds to a stage beneath a dividing wall where the outgoing vapor stream is split (compare Figure 2c). This is why an outgoing vapor stream and an ingoing liquid stream connect this column section to another column section. If stage n s describes a stage above a dividing wall, it would be an additional ingoing vapor and outgoing liquid stream. Finally, if n s corresponds to the column feed stage there would be no additional vapor stream but only one ingoing liquid stream (the column feed). These streams described by the optimization variables can be considered to be known in this section. To receive all process variables, we need to calculate all inner streams and the temperatures within the stages.
After the early work by Lewis and Matheson 1932 and Thiele and Geddes 1933 many methods have been described in the literature to determine all process variables of a conventional distillation column. They include the so-called θ- method Holland 1963;Billingsley 1970;Haas, Alejandro, and Holland 2007, the sum rate method Sujata 1961; Burningham and Otto 1967;Lucia and Li 1992 and the relaxation method Rose, Sweeny, and Schrodt 1958;Ketchum 1979;Mori et al. 1990. Also Newton type methods have been used to solve the corresponding MESH equations Naphtali and Sandholm 1971; Ishii and Otto Figure 4: Theoretical stages of one column section of a mDWC. Increasing control volumes (dashed rectangles). Streams described by optimization variables (thick lines).
1973; Vickery and Taylor 1986. These works aim at solving all MESH equations and determining all process variables. An alternative to using the original MESH equation is to use shortcut methods. Various simplifying assumptions have been developed and repeatedly used in the literature. These include the assumption of constant volatilities, as used in the classic work of Underwood 1948 for a convential distillation column and for dividing wall columns by Amminudin et al. 2001 andHalvorsen andSkogestad 2003a;Halvorsen and Skogestad 2003b, or the assumption of constant molar overflow used in the procedures presented in Levy, van Dongen, and Doherty 1985;van Dongen and Doherty 1985;Julka and Doherty 1990;Zhang and Linninger 2004. These shortcut methods can either be used to approximate the minimum required energy, or their results can be used to calculate initial values for a rigorous optimization In our approach, we aim to solve the original MESH equation directly and make no simplifying assumptions. However, we do not need to solve all MESH equations of the column section as a part of the equations can be included in the optimization problem given in Equation (9) and will be solved by an optimization routine. Moreover, the process variables have to be determined multiple times during the optimization of a mDWC. This means that the equations must be solved reliably and fast. For a conventional distillation column it was previously shown that this can be done by stage-to-stage calculations based on fixed point iterations (see Hoffmann et al. 2017). In the following we will show how stage-to-stage caluclations can be used here to obtain all process variables describing a column section.
As we already know the liquid concentration x 1 , we can use the equilibriumcondition given in Equation (3) and the summation-law to calculate the vapor concentration y 1 and the temperature T 1 in stage 1. We next consider the control volume consisting only of the first stage (depicted in Figure 4 by a dashed rectangle around stage 1). The variables that are still unknown in this control volume are V 1 , L 2 , x 2 . If we start with an initial guess of V 1 we can calculate the liquid concentration using the material-balance in the control volume: where we used L 2 = L 1 +V 1 −V 0 . Using the enthalpy balance of the considered control volume and previously calculated values for x 2 we can calculate a new value for V 1 : Applying Equation (10) and Equation (11) can be repeated until convergence is reached. Using this fixed point iterations we thus obtain values for V 1 , L 2 and x 2 that satisfy both the material-balance and the enthalpy balance for stage 1.
As we obtain a value for x 2 we can continue with the next stage and calculate again y 2 with the use of the equilibrium-condition. We can extend the control volume by one stage (depicted in Figure 4 by a dashed line around stage 1 and 2). This way we receive fixed point iterations for the variables V 2 , L 3 , x 3 . In general for a stage k with 1 ≤ k < n s the fixed point iteration is given for After the fixed point iterations converge to a fixed point s * the missing values can be calculated as follows: If the control volume in Figure 4 is extended to include the stage n s the material-balance and the enthalpy-balance change, as the control volume then includes two further streams. The fixed point iteration can be adjusted accordingly. We will see in the next subsection that an easy condition can be formulated that guarantee that the calculated values in Equation 12 for the next concentration x k+1 are non-negative. This is not the case for the adjusted fixed point iteration. If the values for the liquid concentrations become large negative values the solution of the equilibrium-condition may fail and with it the complete stage-to-stage calculations. To avoid this we stop after stage n s −1 and calculate the remaining streams downwards instead of continuing upwards.
In the downwards calculation we consider a stage index k with n s < k ≤ n. We know already the value of y k . Solving the equilibrium-condition provides values for the liquid concentrations x k . It remains to determine values for y k−1 , V k−1 , L k . Analogously to the above described upwards calculations a fixed point iteration for the downwards calculations determining L k in stage k with n s < k ≤ n is given for s ∈ [0, ∞) by: Again after the fixed point iterations converge to a fixed point s * the missing values can be calculated: The convergence of the fixed point iterations for the upwards and the downwards calculation have been studied for conventional distillation columns in detail by Hoffmann et al. 2017.
While the upwards stage-to-stage calculations solve all MESH equations in the stages 1 to n s − 1, the downwards calculations solve all equations in the stages n s + 1 to n. The equations corresponding to stage n s haven't been considered yet. Also Equation (6) and (7) that guarantee that the vapor or liquid outgoing streams of a split stage have the same concentrations are not yet considered. These equations are included in the optimization problem described at the beginning of this section and are solved by the optimization routine. After convergence of the optimization all MESH equations describing the multiple dividing wall column are fulfilled.

Ensuring positivity constraints
Before we turn to the numerical examples we discuss a property that is important for the stage-to-stage calculations. In every stage of the upwards calculation we compute values for a liquid concentration. To continue with the next stage we have to solve the equilibrium-condition. But so far the fixed point iteration can not guarantee that the calculated values for the liquid concentrations are all positive. If they become large negative values, the equilibrium-condition may not have a solution and we cannot continue with the next stage. This means that the stage-to-stage calculations can fail in the case of negative values for the concentrations. In the remainder of this section we will discuss properties that ensure positiveness of these values.
In the upwards calculation values for the liquid concentrations are given by Equation (12). If we assume that for a stage index k the values for the liquid concentrations x k are non-negative, the values for the vapor concentrations y k calculated using extended Raoult's law are also non-negative. If we now additionally assume that for every component i we also have L 1 ≥ V 0 and the values calculated with Equation (12) are for every s ∈ [0, ∞) non-negative. This means that non-negative values for the concentrations can be guaranteed in every iteration and thus also for the limit.
Inductively the same follows for all k with 1 ≤ k ≤ n s − 1. Analogously, under the assumption that: all values for concentrations calculated in the downwards procedure will be nonnegative.
The inequalities above guarantee that for any choice of the optimization variables all calculated values for the concentrations are non-negative. There is a second possibility to ensure non-negative values. The equations used to derive the fixed point iterations hold for process variables that satisfy all MESH equations, including the equations that were added to the optimization problem. If we start with the corresponding optimization variables, the upwards an downwards calculations will reproduce the process variables. This means that for this optimization variables all values are non-negative. Now, as long as the upwards and downwards calculation are evaluated not too far away from process variables that satisfy all MESH equations, no negative values will occur.
These consideration motivate to use a two step procedure. In a first step both inequalities given in Equations (13) and (14) are added to the optimization problem given in Equation (9). During the optimization the up and downcalculations can always be performed. This means that all functions can be evaluated and the optimization routine can find a solution. As we have added two additional inequalities, the calculated solution is in general not optimal for the original problem without these two additional assumptions. This is why a second step is performed in which we solve the optimization problem without the additional inequalities. As a starting point we use the optimization variables calculated in the first step. These variables already satisfy all MESH equations and the up and down-calculations can be evaluated at this starting point. Whenever a point is reached for which any of the up and down-calculations fails a smaller step can be taken. This way the optimization routine stays close to a feasible point. In our experiments presented in Section 5 the approximation calculated in the first step was always already close to the overall optimal solution.

Details of implementation
To solve the problems described in the next section the new approach has been implemented. As a comparison a black box simulation software and an outer loop for optimization has been implemented.
New Approach: The thermodynamic models, the stage-to-stage calculations and the optimization have been implemented using MATLAB R , version 9.3 2017. More information about the used models and the parameters for the examples can be found in the A. For the models (enthalpy, vapor pressure, activity coefficients) implemented in MATLAB R the parameters are taken from the commercial flowsheet simulater Aspen Plus R , version 11 2019. The models have been validated using example data obtained with Aspen Plus R . The calculations with the new method. have been performed on a 64-Bit Windows computer with 16 GB of RAM and an Intel R Core TM i7-8665U processor. The time was measured using the MATLAB R -function timeit.
To solve the optimization problems we used sequential quadratic programming (sqp) (see for example Nocedal and Wright 2006). More precisely the implementation provided by the MATLAB R Optimization Toolbox TM , version 8.0 and the function fmincon was used. The algorithm was run with the default settings. The accuracy demanded for the final constraint violation was set to 10 −6 . Optimization algorithms can benefit, if clean analytic derivatives are provided and they do not have to be calculated using finite differences. In the stage-to-stage calculations presented in Section 3.1 the MESH equations for a control volume are solved and the dependent variables are calculated. Derivatives for the dependent variables can thus be calculated using the implicit function theorem Königsberger 2002. In MATLAB R the code is compiled after the calculations are started. This can make the calculations slower than necessary. To avoid this loss in speed, the stage-to-stage calculations have been precompiled using MATLAB R Coder TM , version 3.4. This decreases the time needed for each evaluation in the optimization routine.
Aspen Plus R as simulation software: The calculations presented in the following are compared with results obtained using Aspen Plus R , version 11 2019. In these calculations the software Aspen Plus R is considered as a black box simulator solving all MESH equations. The simulation was optimized in an outer loop. This means that in each iteration the optimizer sets all specifications which serve as optimization variables and the simulation is evaluated for these specifications. After a converged simulation the objective function (reboiler duty) and the constraints (purity demands) can be evaluated. Depending on their values the specifications are modified until a satisfying result is found. The procedure is explained in more detail in the works by Kurnatowski et al. 2017 andRänger et al. 2020.
In the present work we used the reboiler duty, the flowrates of the top stream and side product streams as well as the ratio of the splits at the dividing walls as specifications. We again used sequential quadratic programming to solve the optimization problems. Here we worked with two implementations from the Schittkowski optimization suite. Namely the solvers NLPQLP Schittkowski 2009 and MISQP Exler and Schittkowski 2007. The demanded final accuracy was set to 10 −3 . Unfortunately the simulation software Aspen Plus R does not provide derivatives which is why we used derivatives obtained using finite differences. The calculations with Aspen Plus R have been performed on a 64-Bit Windows computer with 16GB of RAM and an Intel R Core TM i5-6500 processor.

Numerical Results
In this section we apply the developed methods to two examples. In a first example we optimize a column with a simple dividing wall splitting a ternary mixture. We will compare the results with Aspen Plus R and show that our implementation is very close to the implementation of Aspen Plus R . The second example is more challenging, since it is about the separation of a quaternary mixture. We demand higher purities and consider two dividing walls. For this example the benefits of the new method become more apparent.

Simple dividing wall column
As a first example we consider the separation of the mixture Benzene, Toluene and p-Xylene. The vapor liquid equilibrium is modeled using the NRTL-model. The corresponding parameters are given in Appendix A. As there are three components this mixture can be separated with one dividing wall. We consider a column with a height of 40 theoretical stages. The dividing wall reaches from stage 11 to stage 30. The feed is located on stage 21, where the stages are again counted from bottom to top starting with 1. We consider a feed flow rate of 3 kmol/h and the concentration of each component is 1 3 mol/mol. We assume a constant pressure of 1 bar within the column. We try to minimize the heat dutyQ R and demand the purities of 0.95 mol/mol in the product streams. The corresponding Petlyuk sequence is shown in Figure 5. The product streams are denoted by A-C. The streams connecting the different sections of the column are numbered by 1 to 10. We will first present the results obtained with the new approach and then compare them to results obtained with Aspen Plus R as black box simulation software.
Originally the problem has around 550 process variables and MESH equations. For a given feed and constant pressure in the column, there remain 5 degrees of freedom. As described in Section 3 the optimization variables need to describe the streams connecting the different sections of the distillation column. For the streams 3, 4 and 7-10 we used the molar flow rates of each component to describe the streams. Instead of describing streams 1, 2 and 5, 6 we considered the molar flow rates of each component of the product streams A and C, the reboiler duty and the condenser duty. From these the streams 1, 2 and 5, 6 can be calculated using the material and enthalpy-balance. Finally, we chose the percentage of product withdrawn at stream B as optimization variable. This way we can also compute the liquid input stream at the top of section C 1,2 from stream 3. Using the new approach the number of variables are reduced to 27  Table 1. An initial heat and condenser duty can be calculated using the enthalpy-balance around the reboiler and the condenser respectively. As pure components are assumed in the product streams, the flow rate for each product stream must be 1 kmol/h. In the optimization routine we used very rough bounds on the optimization variables. The flow rates for each component of a stream are bounded from below by 0 and from above by 100 kmol/h. Finally, the heat duty and the condenser duty were assumed to lie in the intervals [10 kW, 1000 kW] and [-1000 kW, -10 kW] respectively.
As described before in Section 3.2, the optimization was performed in two steps. First, we added the inequalities described by Equation (13) and Equation (14) as constraints to the optimization problem. The algorithm terminated after 41 iterations which took 9.48 seconds. The results are summarized in Table 2. In the table the flow rateṅ and the concentrations of all product streams are given. Furthermore, we give the values for the heat duty and the liquid and vapor split ratios at the dividing wall l s , v s . The value l s denotes the ratio of liquid that flows to the right of the dividing wall on the top of the dividing wall. The value v s denotes the ratio of vapor that goes to the right side of the dividing wall on the bottom of the dividing wall. The results show that the initial values are a good guess for the optimized values. The initial heat duty is approximately 4% higher than the final heat duty. The vapor and liquid split are smaller than predicted. Which means that in an optimal solution a smaller percentage of liquid and vapor flows to the right of the dividing wall. The smallest flow rate of a product stream is attained for the product stream B.
After this first run we dropped the additional inequalities and optimized again starting from the solution found so far. The algorithm terminated after 16 additional iterations which took 4.27 seconds. The results are also shown in Table 2. As removing an additional inequality increases the size of the search space, a better solution is found in the second step. One can see in the table that the flow rate at the top stream becomes a little bit smaller and the purity of the head-product increases. The split ratios differ slightly, the difference in the heat duty is almost negligible. This means that both solutions differ not too much and the solution found in the first step is a good starting point for the second step in which the full search space is explored.
The performed two step procedure was motivated in Section 3.2 to ensure that the stage-to-stage calculations can always be performed. If the optimization problem is directly solved without additional iequalities, the stage-to-stage calculations may be unstable. If they fail too often, the solver may not converge to an optimal or even feasible solution. In this easy example the initial values calculated with the V min diagram are already quite close to the final solution. This is why we tried to solve the optimization problem directly without the additional inequalities. The solver had some difficulties but converged after 102 iterations (65.25 seconds) to the same optimal solution. In 29 of the iterations the stage-to-stage calculations failed which was handled by a forced step back in the line-search of the sqp-algorithm. In this example, it is in principle possible to avoid the two step procedure. However, the solver needed more iterations and more time. The failures in the evaluations show that it is more difficult to solve the problem directly and it can beneficial to follow the proposed two step procedure.
We next compare the results of the new method with results obtained with Aspen Plus R and an optimization in an external loop. As described before there are 5 degrees of freedom left, which can be chosen as optimization variables. We chose the flow rates of the product streams A and B (both ranging from 0.5 to 1.5 kmol/h), the reboiler duty (25 to 200 kW) and the two split ratios l s (0.2 to 0.9) and v s (0.15 to 0.9) as optimization variables. The objective was again to minimize the reboiler dutyQ R and we demanded the same purities of 0.95 mol/mol in each product stream. The optimization routine needed 18 iterations which took 534.74 seconds. This is about 39 times more time than the new approach needed. Both results, the one obtained with the new approach and with Aspen Plus R as simulation software, are shown in Table 3. Again we list the product streams, the heat duties and the split ratios. Comparing both solutions several observations can be made. First both solutions are close together. The reboiler duty calculated with the new approach is smaller, but only by 0.29%. The difference in the reboiler duties is larger than the difference between both solutions calculated with and without additional inequalities (see Table 2). The purity demands are slightly violated in the solution with the Aspen Plus R simulation. The largest violation in the solution calculated with the new approach is in the order of 10 −16 . In the new approach the purity demands can be calculated as a linear constraint on the optimization variables. These linear constraints are respected by the sqp-algorithm up to a numerical precision. The situation is different using a simulation software together with an external optimizer. The concentrations are outputs of the simulation and their dependency to the optimization variables is highly nonlinear. The largest difference in the solutions can be observed for the split ratios. This means that in this example the solution is not very sensitive to changes in the split ratios, which is known from literature Halvorsen and Skogestad 2001. There are two possible reasons for the different results found by both approaches. First, the optimization performed with an outer loop terminates earlier with a slightly sub-optimal solution. Second, the implemented models differ. If the models differ by too much the results are not comparable. To exclude the second possibility we determined the difference between a simulation in Aspen Plus R and our own implementation. After the new approach converged we can calculate with the final optimization variables all process variables by evaluating the map given in Equation (8). This includes values for the reboiler duty, the split ratios and flowrates of the product streams. These values can then be used as specifications in a simulation performed with Aspen Plus R . This way we obtain process variables for the same specification once calculated with the new approach and once calculated using Aspen Plus R . Comparing the different set of process variables show that they are very close together. The largest difference for one of the product concentrations is 4.69 · 10 −6 . This means that the difference in the implementations is negligible. The solution found with the new approach has indeed a smaller violation of the purity specification and has a better objective value than the solution found by an outer optimization.
One of the reasons for the slightly worse solution is the lack of clean analytic derivatives. Close to the optimal solution the derivatives calculated numerically suggest the solver, that he cannot come closer to the solution anymore. This first example already shows that it is possible to optimize dividing wall columns with less iterations than using the simulation software Aspen Plus R and an optimization in an outer loop. Because of the availability of analytic derivatives numerically cleaner solutions can be calculated.

Multiple dividing wall column
We next consider an example that is more difficult to solve. We want to separate a quaternary mixture consisting of Ethanol, n-Propanol, i-Butanol and n-Butanol in a dividing wall column with two dividing walls as shown in Figure 1c. This system was suggested in literature to be suited for the considered mDWC Preißinger, Ränger, and Grützner 2019. Again the vapor liquid equilibrium is modeled using the NRTL-model and the corresponding parameters are given in A. The distillation column has a total height of 78 theoretical stages and is operated at a constant pressure of 1 bar. The separation is performed in a distillation column with two dividing walls. The first dividing wall splits the stages 14 to 52 an the second splits the stages 40 to 65. The second dividing wall is located to the right of the first dividing wall. The feed is located on stage 40. We consider a feed of 100 mol/h and molar concentrations of 1/4 mol/mol for each component. We try to minimize the reboiler dutyQ R and demand the purities of 0.98 mol/mol in the product streams. The different sections of the column organized as Petlyuk sequence are shown in Figure 5. The product streams are denoted by A-D. The streams connecting the different sections are numbered by 1 to 16. As in the previous example, we will first present the results obtained with the new approach and then compare them to results obtained with Aspen Plus R as black box simulation software.
In the presented approach the optimization variables describe the streams connecting the different sections of the distillation column. Analogous to the previous example we chose the following process variables: the flow rates of each component in the streams 3-6 and 9-16; the flow rates in the product streams A and D; the reboiler and the condenser duty and the percentage withdrawn in the side streams B and C. In total with the new approach we obtain an optimization problem with 60 optimization variables and 52 equality constraints that remain to be solved. As an comparison the original simulation problem has around 1600 process variables and MESH equations.
Again starting values for the optimization variables are needed. A possibility would be to use again streams obtained from a V min diagram. In our numerical experiments we were able to find solutions using the new approach starting with these streams. Unfortunately, the optimization using Aspen Plus R as simulator did not converge. The V min diagrams assume an infinite number of stages. In the current example the number of stages is too small so that the approximation made with an infinte number is not very close. The initial simulation performed with values of the V min diagram predicted purities of around 60% in some of the product streams. This is why, it was suggested to adjust the initial flowrates of the internal streams Ränger and Grützner 2021. The flow rates are sufficient to calculate the initial specifications used in the simulation with Aspen Plus R . However, initial values for the concentrations are needed for the new approach. The new approach is not too sensitive to the choice of initial values for the Figure 6: Petlyuk sequence corresponding to a distillation column with two dividing walls. optimization variables, so we did not adjust the concentrations and simply kept the concentrations predicted by the V min diagram. The initial values for the connecting streams used for the results in the following are given in Table 4. As pure components are assumed in the product streams, the flow rate for each product stream must be 25 mol/h. The flow rates for each component of a stream are again bounded from below by 0 and from a above by 10 kmol/h. The heat duty and the condenser duty are assumed to lie in the intervals [1 kW, 10 kW] and [-10 kW, -1 kW].
Using the improved starting values, we can proceed with the actual optimization by introducing the minimization of the reboiler duty as an objective. As before we first used the additional inequalites given in Equations (13) and (14). The calculation terminated successfully after 56 iterations (25.23 seconds). The product streams of the obtained solution are given in Table 5. Additionally, values for the reboiler duty as well as the split ratios are given. Where v s,1 and v s,2 are the ratios of vapor that flows to the right of the first and the second dividing wall. The ratios of liquid that flows to the right of the dividing walls is denoted by l s,1 and l s,2 .
To explore the full search space we removed the additional inequalities and ran the optimization again using the previously calculated values as starting point. The optimization needed 12 iterations (5.31 seconds). The results are given in Table 5. In contrast to the first example studied in Subsection 5.1 the solution with and without additional inequalites differ more. The objective value in the second run is around 7.08% smaller. The difference in the product streams is very small. The results show that in an optimal solution less is withdrawn in the side product streams than in the bottom and distillate stream. We again compare the results of the new approach to results obtained with an external optimization routine and Aspen Plus R as simulation software. As specifications that are varied by the external optimizer we used the reboiler dutẏ Q R (4.4 to 5 kW), the total molar flow rates of the product streams A, B, C (24.5 to 25.5 mol/h) and the four split ratios v s,1 , v s,2 , l s,1 , l s,2 . All of them are limited by 0.1 from below and 0.95 from above, except for v s,2 the lower bound was set to 0.2. The external optimization routine needed 189 iterations which took 11.77 hours. Both results, the one obtained with the new approach and with Aspen Plus R as simulation software, are shown in Table 5. Again we list the product streams, the heat duties and the split ratios. Both solutions are similar. The split ratios only differ slightly this time. The maximal violation of the purity specification for the solution calculated with Aspen Plus R is 0.01%. The new approach found a solution with a reboiler duty that is 1.73% smaller.
In the above calculations we used initial values calculated with V min diagrams and adapted as described in Ränger and Grützner 2021. These values give a good initial estimate and, as shown in Table 5, are already close to the final optimal solution. The main reason for the need of good initial values was the comparison with the calculations performed with Aspen Plus R . To get 5.09 · 10 −20 5.16 · 10 −20 1.87 · 10 −20 x P r /(mol/mol) 0 3.78 · 10 −6 4.13 · 10 −6 3.29 · 10 −6 x Is /(mol/mol) 0 0.0200 0.0200 0.0200 x Bu /(mol/mol) 1 0.9800 0.9800 0.9800 good values with Aspen Plus R , good initial values are needed. In the rest of this section, we want to investigate the impact of these initial values on the convergence of the new approach. Therefore, instead of initial values computed with V min diagrams we consider random initial values constructed as follows: • The reboiler and condenser duty are chosen uniformly distributed in the interval [1, 10 ] kW.
• The molar flow rate of the product streams A and D are chosen uniformly distributed in the interval [0,100] mol/h. The concentrations for these two streams are chosen as in Table 4.
• The molar flow rate of each component in the internal streams 3-6 and 9-16 which have a concentration larger than 0 in Table 4 are chosen uniformly distributed in the interval [0,250] mol/h.
The bounds were chosen in such a way that the intervals contain the values computed with the V min diagrams with and without the modifications. The random initial points are reasonable in the sense that concentrations are set to 0 which are expected to be close to 0 in a feasible solution, but otherwise the random initial values can vary very roughly and can be far away from a solution. With these new initial values we followed the two step procedure by first adding the inequalities given in Equations (13) and (14) and then solving the full optimization problem. We repeated this for 100 randomly chosen initial values. The results are summarized in Table 6. For each run, we checked whether the optimality conditions were met for the final values. The number of these runs divided by the total number of runs performed gives the success rate. We also collected the final objective value and the time needed for each successful run. As a comparison we used ten times the starting values proposed by the V min diagrams (orig. V min ) and the values obtained by the modification introduced by Ränger and Grützner 2021 (mod. V min ). As expected, due to the high non-linearity of the MESH equation, the calculations do not converge from all possible choices of initial values. However, the results show that even if the starting values are chosen randomly, convergence can be achieved for a majority of the test runs. The success rate with random initial values can be further increased with a multi start strategy. For example, the probability of not finding the optimal solution decreases to a value less than 0.01% for 3 independent repetitions. With all three different strategies, we found the same optimal solution. In fact, all successful runs converged to the same solution. This means that although the problem formulation is highly non-linear, no different local solutions were found. Despite the fact that convergence can often be achieved with random initial values, the use of good starting values is still beneficial. First, in our experiments, the good initial values always converged. Second, fewer iterations were needed, resulting in less computation time.
In the calculations with the multiple dividing wall column the benefits of the new approach become evident. First, the stability of the calculations is improved. To obtain the results using Aspen Plus R quite some effort had to be taken. The initial values had to be adjusted, bounds for the specifications had to be chosen carefully in order to avoid non converging simulations and to find an optimal solution. The stage-to-stage calculation used in the new approach converge on a large range of variables. This is why rough bounds could be used. We could show that even with random initial values convergence can be achieved in the majority of cases. Second, with the new approach solutions of a higher numerical quality could be found. The found solution shows a smaller violation of the purity specification and has a smaller objective value than the solution calculated with Aspen Plus R and an external optimization. Finally, the reduction in the calculation times is significant. In total the calculations with the new approach only took 30.55 seconds while the calculation with an external optimizer took more than 11 hours. This means that the new calculations were more than 1000 times faster.

Conclusion
Instead of solving the MESH equations repeatedly and performing the optimization in an outer loop, the presented approach embeds the process simulation into an optimization problem. In this way, the MESH equations are solved only once and the optimization is performed simultaneously. To avoid a large and very non linear optimization problem, not all MESH equations are considered as constraints and not all process variables are considered as optimization variables. A subset of the process variables is chosen as optimization variables. From this subset, the remaining process variables are calculated by a robust and fast calculation procedure. This procedure guarantees that a large subset of the MESH equations is already satisfied.
To reduce the number of variables for multiple dividing wall columns we introduced robust stage-to-stage calculations for sections of the distillation column. In each section all process variables are calculated from top and bottom towards the middle until the splitting or feed stage is reached. In each stage the MESH equations are solved using fixed-point iterations.
In two examples we were able to present the benefits of the new method. In both examples the new method converged easily towards an optimal solution. For the second example we demonstrated the robustness of the proposed approach by using random initial values. The optimality could be confirmed by a calculation performed by the commercial flowsheet simulator Aspen Plus R in combination with an external optimization routine. Comparing both solutions showed that the new approach found a solution of a higher quality and more stable. In terms of the needed time to find an optimal solution, the new approach totally outperformed the classical approach. For a distillation column with two dividing walls the new approach was 1000 times faster.
The developed approach is not limited to multiple dividing wall columns, but can in principle also be applied to other operating units and flowsheets consisting of several operating units. The determination of suitable optimization variables and the reduction of the descriptive system of equations must be done individually for each different type of operating unit. The approach is particularly powerful and outperforms the state of the art by far when there is a large number of process variables and a strongly nonlinear system of equations present.
The gained stability and the speed improvement will allow further investigation in the design of multiple dividing wall columns in future work. With the new approach one does not need to carefully set bounds to avoid non convergent areas. The presented method converges on a large range of variables. Instead of waiting almost half a day for results the time needed for solving one optimization problem, the time needed with the new approach is reduced to around one minute. This means that much more calculations and different optimization problems can be solved in a shorter time.
The used parameters for the different components are given in Table 7 and  Table 8. The parameters are taken from Aspen Plus R , version 11 2019.  For the liquid enthalpy of a mixture with molar concentrations x the mixing enthalpy is added, which is given by: where R is the gas constant and γ i the activity coefficient of the i-th component.

A.2 Vapor pressure
The vapor pressure has been calculated using extended Antoines equation: The units are in Pa. The parameters are given in Table 9.

A.3 Activity coefficients
The ctivity coefficients have been calculated using the NRTL-model Renon and Prausnitz 1968 with τ ij = a ij + b ij /T and the non-randomness parameters α ij . The parameters are taken from Aspen Plus R , version 11 2019. For the first mixture they are given in Table 10. For the second mixture they are given in Table 11 Table 10: NRTL-model parameters for the description of the vapor-liquid equilibrium in the ternary system Benzene, Toluene and p-Xylene.