**A Hybrid Multi-Objective Optimization Framework for Preliminary Process Design Based on Health, Safety and Environmental Impact**

**Shin Yee Teh 1, Kian Boon Chua 1, Boon Hooi Hong 1, Alex J. W. Ling 1, Viknesh Andiappan 2, Dominic C. Y. Foo 1, Mimi H. Hassim <sup>3</sup> and Denny K. S. Ng 1,2,\***


Received: 25 February 2019; Accepted: 25 March 2019; Published: 8 April 2019

**Abstract:** Due to increasingly stringent legal requirements and escalating environmental control costs, chemical industries have paid close attention to sustainable development without compromising their economic performance. Thus, chemical industries are in need of systematic tools to conduct sustainability assessments of their process/plant design. In order to avoid making costly retrofits at later stages, assessments during the preliminary design stage should be performed. In this paper, a systematic framework is presented for chemical processes at the preliminary design stage. Gross profit, Health Quotient Index (HQI), Inherent Safety Index (ISI) and the Waste Reduction (WAR) algorithm are used to assess the economic performance, health, safety and environmental impact of the process, respectively. The fuzzy optimization approach is used to analyse the trade-off among the four aspects simultaneously, as they often conflict with each other. Deviation between the solution obtained from mathematical optimization model and process simulator is determined to ensure the validity of the model. To demonstrate the proposed framework, a case study on 1, 4-butanediol production is presented.

**Keywords:** input-output model; fuzzy optimization; process synthesis; preliminary stage design

#### **1. Introduction**

Process design is a core element in the field of chemical engineering. It can be considered a centre point, bringing together all chemical engineering components as a whole. Process design is associated with creating processes or improving existing processes. An integral part of process design is *process synthesis*. Process synthesis is defined as "the discrete decision-making activities of conjecturing (1) which of the many available component parts one should use, and (2) how they should be interconnected to structure the optimal solution to a given design problem." [1]. The field of process synthesis has seen significant developments since its inception in the 1960s led by the late Roger W. H. Sargent [2]. Most notably, key pioneering contributions have been an integral part of establishing what process synthesis is and what it entails. For instance, Nishida et al. [3] provided an important overview of developments within the boundary of process synthesis and defined process synthesis as "an act of determining the optimal interconnection of processing units as well as the optimal type and design of the units within a process system". Meanwhile, the *Onion model* [4] was proposed as a

systematic overview and guide to process synthesis thinking. The Onion model emphasized that a reactor is designed first followed by separation and recycle streams, heat recovery systems and utility systems [4]. Foo and Ng [5] then extended the *Onion model* further by incorporating material recovery and treatment systems (see Figure 1). On the other hand, Douglas' model is another well-accepted decision hierarchy approach that was proposed for process synthesis in the late 80s [6].

**Figure 1.** The extended Onion Model [5]. Reproduced with permission from Foo, D.C.Y.; Ng, D.K.S., Process Integration for Cleaner Process Design; In Handbook of Process Integration: Minimisation of Energy and Water Use, Waste and Emissions; published by Elsevier Science, 2013.

To date, there has been a vast number of developed process synthesis methodologies. Often, these methodologies pay close attention to the technical and economic performance of a process design. Meanwhile, other aspects such as safety, occupational health and environmental impacts are typically considered at the later/final stages of design [7,8]. However, several papers argue that these aspects should be considered at the preliminary stage of process synthesis as the cost of process improvement and operational risks can be significantly lowered compared to at the later stages (detailed design) [9,10]. In this respect, several alternative methods for evaluating occupational health, safety and environmental impacts have been presented.

In the area of safety evaluations, Edwards and Lawrence [11] developed the earliest method for assessing inherent safety in a given process. This method is known as *Prototype Index of Inherent Safety* (PIIS). PIIS ranks the inherent safety level of alternative chemical process routes based on main reactions and parameters such as pressure, temperature, yield, heat of reaction, inventory, flammability, toxicity and explosiveness. Despite considering the mentioned parameters, PIIS focuses solely on the main reaction of a process and not the other parts of the process. In view of this, Heikkila [12] proposed an alternative approach called the *Inherent Safety Index* (ISI). ISI considers the same factors as in the PIIS along with additional ones such as corrosions, side reactions, and inventory for both inside and outside the battery limits, types of equipment and process structure. Later, Palaniappan et al. [13]

added a few supplementary indexes such as worst reaction index, overall chemical index and total chemical index to the previously proposed ISI to form *i*-Safe.

Aside from safety, various methods have been developed to evaluate occupational health. For example, Hassim and Edwards [14] developed Process Route Healthiness Index (PRHI) to rank the process alternatives based on the potential of working activities and process conditions that may harm workers. However, PRHI is not suitable for assessment in the preliminary design stage, as it requires complete process information (e.g., points for manual handling etc.). In this respect, Hassim et al. developed the Inherent Occupational Health Index (IOHI) [10], which is suitable for application at the initial stage of the process research and development (R&D). This work was then extended to a more detailed assessment called the Health Quotient Index (HQI), which uses detailed process information available in the preliminary design stage [15]. HQI is able to rank the process alternatives based on health risk values from the fugitive emissions and to calculate the risk of a given process.

Apart from occupational health, there are many well-established approaches reported in the literature on environmental impact assessments. The Environmental Hazard Index (EHI) [16] and Atmospheric Hazard Index (AHI) [17] are among the early methods that assess inherent environmental performance of chemical process routes. Later, Gunasekera and Edwards [18] integrated EHI with AHI into a new method which is known as Inherent Environmental Toxicity Hazard (IETH). IETH can estimate the inherent environmental friendliness of a chemical plant in all media including air, soil and aquatic due to total loss of contaminants. Cabezas et al. [19] then introduced the generalized Waste Reduction (WAR) algorithm. This algorithm determines the total environmental impacts based on the Potential Environmental Index (PEI) balance, which assigned environmental impact values to different pollutants. Young and Cabezas [20] extended PEI balance to integrate the energy consumption of the chemical process into the environmental evaluation. Similarly, Andiappan et al. [21] presented the incremental environmental burden assessment (IEBA), which used the concept of economic potential assessment method developed by Seider et al. [22] to determine the new environmental burden for each process.

It is essential to note that the abovementioned evaluation methods focus only on assessing the performance of a process design based on just one aspect (e.g., safety or occupational health or environmental impact) and not addressing them simultaneously. In this respect, there have been attempts to develop process synthesis approaches that address safety, occupational health and environmental impacts simultaneously. For instance, Azapagic et al. [8] and Othman et al. [23] presented sustainable assessment and design selection approaches that consider economic, environmental, and social aspects. Besides, Al-Sharrah et al. [24] presented a multi-objective optimization model which considers environmental impact, economic performance and operational risk simultaneously for the petrochemical industry. Liew et al. [25] developed a systematic approach to the screening of sustainable chemical reaction pathways at the research and development (R&D) stage. Based on the proposed approach, fuzzy optimization is used to trade-off economic performance, health, safety and environmental impact. Similarly, Ng et al. [26] presented a multi-objective process synthesis approach which considered economic performance, health, safety and environmental impact for biorefineries. Following this, a visualization tool called the *Piper diagram* was also proposed for considering economic, safety and environmental aspects simultaneously [27]. However, due to its graphical nature, the Piper diagram is unable to consider process optimization.

It is clear that the abovementioned work presents approaches to evaluate safety, occupational health and environmental impacts at the preliminary stage of design. However, these approaches focus on preliminary screening without considering operating conditions in a given process. The choice of process operating conditions has a process direct impact on screening decisions and this should be given closer attention. To address this, process simulators can be used to analyse operating conditions for a given process. Process simulators can be used to represent chemical processes in terms of mathematical models and solve them to attain insights on their performance [28]. Process simulation is an essential and complementary task of process synthesis, as it predicts how a process design would behave under defined operating conditions [29,30]. Despite the benefits of process simulators, they suffer from various limitations. For instance, simulators are limited to considering a single objective (e.g., product yield or production rate) at a given run. This limitation disables simulators from finding a trade-off between multiple aspects such as process economics, occupational health, safety and environmental impacts for a given process design. In order to consider multiple aspects simultaneously, multi-objective optimization is required.

As such, this paper presents a systematic framework which combines the benefits from both process simulation and multi-objective optimization to address economic, environmental, health and safety aspects simultaneously at the preliminary design stage. This paper is structured as follows: First, a formal problem statement is given in Section 2. A systematic framework for preliminary process design is presented in Section 3; To illustrate the proposed framework, a case study on the screening of 1,4-butanediol production processes is presented in Section 4; and the work is finally concluded in Section 5.

#### **2. Problem Statement**

The problem addressed in this work is stated as follows; given a set of alternative processes *k* to produce chemical product (in a given output stream) *p*. Each alternative process *k* has a set of unit operation *j* with operating capacities, *xj,k* and a set of (output) stream *p* to or from unit operation *j* represented by matrix a*p,j,k*. Alongside this, each alternative process *k* differs in economic performance, health, safety and environment impacts. In this work, four objectives are considered for performing a preliminary evaluation and screening of alternative process *k*. The economic performance is determined based on the gross profit of alternative process *k* (GP). The health impacts of process alternative *k* are evaluated using Health Quotient Index (HQI). Meanwhile, the safety and environmental impact of each alternative process *k* are determined via Inherent Safety Index (ISI) and Potential Environmental Index (WAR) respectively.

This work combines the advantages of process simulation, input-output modelling (IOM) and fuzzy optimization to determine the trade-off between the GP, HQI, ISI and WAR objectives. In fuzzy optimization, the degree of satisfaction (λ*k*) is introduced to quantify the degree of satisfaction for four objectives in each alternative process *k*. Following this, the alternative process *k* with the highest λ*<sup>k</sup>* will be selected. A detailed description of the proposed framework to address the stated problem is presented in the following section.

#### **3. Framework for Preliminary Process Design**

Figure 2 presents the proposed framework for preliminary process design, which considers economic performance, health, safety and environment impacts simultaneously. As shown, the proposed framework begins by identifying the product *i* that is expected to be produced. Next, alternative process *k* which produce product *i* is determined. In order to analyse the performance of each alternative process *k*, process simulation tools (e.g., Aspen HYSYS, SuperPro Designer, PRO/II, etc.) can be used to simulate the process. Based on the simulation results, process data such as mass flow rates and energy requirements for each process unit *j* can be determined. These data are then used to develop an input-output model for each process. Input-output model was presented by Leontief [31,32] to analyse the relationship among the raw materials requirement, goods production and the exchange of materials within different economic sectors. Tan et al. [33] extended the usage of input-output model in life cycle assessment and ecological footprint analysis. Aviso et al. [34] integrated the approach with optimization framework to analyse the eco-industrial supply chain under a water footprint constraint. Recently, the input-output model has been used to describe the material and energy balance of processes in a system [35]. Based on the approach [35], the optimal operational adjustment in multi-functional energy systems in response to process inoperability can be determined. Kasivisvanathan et al. [36] also used the input-output model with robust optimization for process synthesis and design of multi-functional energy systems with uncertainties. Most recently,

Foong et al. [37] extended the use of the input-output model for sustainable oil palm plantation development. A detailed description on the procedure to develop an input-output model is presented with an example in Andiappan et al. [38]. The input-output model for each alternative process *k*, given as in Equation (1)

**Figure 2.** Systematic design framework for preliminary process design.

In Equation (1), a*p,j,k* represents the matrix of input or output mass flow rate of stream *p* to/from unit operation *j* in process *k*. *xj,k* represents the capacity or size of the unit operation *j* in process *k* and *yp,k* is the net output of stream *p* in process *k*. A positive value of *yp,k* indicates that it is an output stream which is either product or effluent stream whereas a negative value indicates that it is an input stream and a value of zero denotes that it is an intermediate stream.

Meanwhile, the lower and upper limits of the capacity of unit operation *j* in process *k* are presented in Equation (2):

$$\mathbf{x}\_{j,k}^{\mathcal{L}} \le \mathbf{x}\_{j,k} \le \mathbf{x}\_{j,k}^{\mathcal{U}} \forall j \forall k \tag{2}$$

In addition to Equations (1) and (2), the economic, environmental, health and safety impact for each process *k* is determined. Firstly, in this work, economic performance for each process *k* is measured via gross profit as shown in Equation (3):

$$GP\_k = \left(\sum\_{p=1}^P y\_{p,k}^{Prod} \mathbb{C}\_{p,k}^{Prod} + \sum\_{p=1}^P y\_{p,k}^{Raw} \mathbb{C}\_{p,k}^{Raw} \right) \text{AOT}\_k \forall k \tag{3}$$

where *y*Prod *<sup>p</sup>*,*<sup>k</sup>* and *<sup>y</sup>*Raw *<sup>p</sup>*,*<sup>k</sup>* are the mass flow rate of the products and raw materials in process *k* respectively; CProd *<sup>p</sup>*,*<sup>k</sup>* and CRaw *<sup>p</sup>*,*<sup>k</sup>* are the product price and cost of raw materials in process *k* respectively; AOT*<sup>k</sup>* is the annual operating time for process *k*. Note that *y*Raw *<sup>p</sup>*,*<sup>k</sup>* is obtained from the input-output model as a negative value. This is because it is an input stream into the system, as stated previously in Equation (1).

On the other hand, the environmental impact of each process *k* is also evaluated. For this work, the WAR algorithm is adapted to evaluate the environmental impact because it can determine the average possible impact of a chemical process on the environment, based on environmental impact values of different pollutants and their respective mass flows. According to the WAR algorithm, the total environmental impact generated by process *k*, *WARk* can be expressed as in Equation (4) as follows.

$$\mathsf{WAR}\_{k} = i\_{k}^{(\mathsf{cp}),\text{ out}} + i\_{k}^{(\mathsf{cp}),\text{ in}} + i\_{k}^{(\mathsf{cp}),\text{ out}} \mathsf{W}k \tag{4}$$

where *i* (cp),in *<sup>k</sup>* and *i* (cp),out *<sup>k</sup>* are the input and output rates of impact in the process *k*, *i* (ep),out *<sup>k</sup>* is the power consumption for process *k*. The value of *i* (cp),out *<sup>k</sup>* , *i* (cp),in *<sup>k</sup>* and *i* (ep),out *<sup>k</sup>* are determined through Equations (5)–(9) respectively:

$$\mathbf{y}\_k^{(\text{cp}),\text{out}} = \sum\_{i=1}^{I} y\_{i,k}^{Out} \sum\_{a=1}^{A} \sum\_{j=1}^{I} w\_{a,j,k} \text{PE}I\_a \mathsf{V} k \tag{5}$$

$$\mathbf{u}\_k^{(\text{cp}),in} = \sum\_{i=1}^I y\_{i,k}^{I\text{n}} \sum\_{a=1}^A \sum\_{j=1}^J w\_{a,j,k} PEI\_a \forall k \tag{6}$$

$$\mathbf{y}\_k^{(\text{ep}),\text{out}} = \sum\_{i=1}^{I} -y\_{i,k}^{\text{Elec}} P E I^{\text{Elec}} \mathbb{V} k \tag{7}$$

$$PEI\_a = \sum\_{l=1}^{L} \alpha\_l \text{PEI}\_{a,l} \forall k \tag{8}$$

$$PEI^{\text{Elec}} = \sum\_{l=1}^{L} \alpha\_l \text{PEI}\_l^{\text{Elec}} \forall k \tag{9}$$

where *y*Elec *<sup>i</sup>*,*<sup>k</sup>* is the power consumption of process *<sup>k</sup>*; *PEIa* and *PEI*Elec are the score of the potential environmental impact of chemical component *a* and electricity respectively; α*<sup>l</sup>* is the weighting factor of the impact at category *l*; PEI*a,l* and *PEI*Elec *<sup>l</sup>* are the potential environmental impact score of chemical component *a* and electricity at each of the category *l* respectively. Note that a total of eight categories are considered for PEI, that is, human toxicity potential by exposure, both dermal and inhalation (HTPE), human toxicity potential by ingestion (HTPI), aquatic toxicity potential (ATP), terrestrial toxicity potential (TTP), ozone depletion potential (ODP), global warming potential (GWP), acidification potential (AP) and photochemical oxidation potential (PCOP).

In order to evaluate the health impact of process *k*, the Hazard Quotient Index (HQI) is selected for this work. This is because HQI can be used for simple process flow diagrams (PFDs) containing limited information, namely process drawings and process descriptions. Such method is suitable for the case of preliminary process design optimization as it allows for the comparison of alternative processes by ranking them based on the risk value based on minimal available information. The calculation of HQI consists of four parts, i.e., estimation of fugitive emissions, air volumetric flow rate, airborne chemical concentration and the health risk (HQI) value. Fugitive emission of each of the chemical components *a* in process *k*, *ma,k* can be calculated via Equation (10):

$$m\_{a,k} = \sum\_{j=1}^{l} \sum\_{p=1}^{P} \mathbf{x}\_{j,k}^{0.5} \mathbf{w}\_{a,j,k} \text{FE}\_{p,j,k} \forall a \forall k \tag{10}$$

where *xj,k*, w*a,j,k* and FE*p,j,k* are the capacity of the unit operations *j*, weight composition of chemical component *a* and the estimated fugitive emission rates in stream *p* to or from unit operation *j*respectively. The pre-calculated fugitive emission rates database for the unit operation stream can be obtained from Hassim et al. [15]. The air volumetric flow rate is determined by the following Equations (11)–(14).

$$\mathbf{A}\_k^T = \sum\_{j=1}^J \mathbf{A}\_{j,k} \forall k \tag{11}$$

$$\mathbf{s}\_k = \left(\mathbf{A}\_k^T\right)^{\frac{1}{2}} \forall k \tag{12}$$

$$\mathbf{A}\_k^c = h\_k \mathbf{s}\_k \forall k \tag{13}$$

$$\mathbf{Q}\_{k} = \mathbf{v} \mathbf{A}\_{k}^{c} \mathbf{V} k \tag{14}$$

where A<sup>T</sup> *<sup>k</sup>* and A*j,k* are the total process floor area and the floor area of each unit operations *j* in process *k* respectively; *sk* is the side length of the process *k*; Ac *<sup>k</sup>* is the cross-sectional area of the process; *hk* is the average height of the main unit operations in process *k*; Q*<sup>k</sup>* is the air volumetric flow rate; v is the wind speed. The average concentration of the chemical components in the air at the downwind edge of the plot area, *Ca,k* can be determined using Equation (15):

$$\mathbb{C}\_{a,k} = \frac{m\_{a,k}}{\mathbb{Q}\_k} \forall a \forall k \tag{15}$$

HQI of component *a* in process *k* can be calculated using Equations (16) and (17):

$$HQI\_{a,k} = \frac{\mathbb{C}\_{a,k}}{\mathbb{C}\_{a}^{\text{EL}}} \forall a \forall k \tag{16}$$

$$HQI\_k = \sum\_{a=1}^{A} HQI\_{a,k} \forall k \tag{17}$$

where *HQIa,k* is the HQI of each chemical components *a*; CEL *<sup>a</sup>* is a constant which represents the threshold limit of the chemical component *a*; *HQIk* is the total HQI of process *k*.

The safety impact of a process is assessed via Inherent Safety Index (ISI). ISI is selected because it is a simple scoring method that can be incorporated into an optimization framework. Moreover, ISI considers all parts of the process and equipment, unlike the PIIS method. The total ISI score for

process *k* is represented by *ISIk*. *ISIk* is calculated by summing the Chemical ISI, *Ik* CI and Process ISI, *Ik* PI as shown in Equation (18):

$$ISI\_k = I\_k^{\rm CI} + I\_k^{\rm PI} \forall k \tag{18}$$

The sub-index for Chemical ISI, *Ik* CI is expressed as below:

$$I\_k^{\text{CI}} = I\_k^{\text{RM, max}} + I\_k^{\text{RS, max}} + I\_k^{\text{INT, max}} + I\_k^{\text{FL, max}} + I\_k^{\text{EX, max}} + I\_k^{\text{TOX, max}} + I\_k^{\text{COR, max}} \forall k \tag{19}$$

where *I* RM, max *<sup>k</sup>* , *I* RS, max *<sup>k</sup>* , *I* INT, max *<sup>k</sup>* , *I* FL, max *<sup>k</sup>* , *I* EX, max *<sup>k</sup>* , *I* TOX, max *<sup>k</sup>* and *I* COR, max *<sup>k</sup>* are the sub-index for the factor of the heat of main reaction, heat of side reaction, chemical interaction, flammability, explosiveness, toxic exposure and corrosiveness of the chemical components present in the process *k* respectively. The sub-index for process ISI, *Ik* PI is given as:

$$I\_k^{\rm PI} = I\_k^I + I\_k^{\rm T, \max} + I\_k^{\rm P, \max} + I\_k^{\rm EQ, \max} + I\_k^{\rm ST, \max} \forall k \tag{20}$$

where *I* I *k* , *I* T, max *<sup>k</sup>* , *I* P, max *<sup>k</sup>* , *I* EQ, max *<sup>k</sup>* and *I* ST, max *<sup>k</sup>* are the sub-index for the factor of inventory, process temperature, process pressure, equipment safety and safe process structure, respectively. The calculations for both *Ik* CI and *Ik* PI are performed on the basis of the worst-case scenario. Among all chemical components present in the process, the scores of the chemical with the most severe hazard in terms of flammability, explosiveness and toxicity are used in the *Ik* CI calculation. Besides, the highest temperatures of the main and side reaction and the worst possible chemical substance interaction in the process are considered. Meanwhile, the maximum expected values for inventory, process temperature, and pressure and the worst process structure are taken into account for the calculation of *Ik* PI. The sub-index for inventory is influenced by the total output of the process, *y*Out *<sup>k</sup>* which is calculated using Equation (21):

$$y\_k^{\text{Out}} = \sum\_{p=1}^{P} y\_{p,k}^{\text{Out}} \forall k \tag{21}$$

where *y*Out *<sup>p</sup>*,*<sup>k</sup>* is the mass flow rate of the output stream for product and effluent streams. The score of the sub-index for inventory is assigned using Equations (22) and (23) [39]. In Equation (22), Lk,*<sup>r</sup>* and Uk,*<sup>r</sup>* are the lower and upper bounds of a given criteria *r*. The criteria *r* can often be in the form of a range. For instance, *r* = 1 may refer to the criteria where mass flowrate is between the range of 300 – 500 kg/hr. In this context, 300 kg/hr would be Lk,*<sup>r</sup>* and 300 kg/hr would be Uk,*r*. If *y*Out *<sup>p</sup>*,*<sup>k</sup>* falls between this range *r* = 1 or meets this criteria, Equation (22) will assign a specific score S*k,r*, which is activated by binary variable *bk,r*.

$$\mathbf{u}\left(\mathbf{L}\_{k,r} - \mathbf{U}\_{k,r}\right) \times b\_{k,r} < y\_{k,r}^{\text{Out}} - \mathbf{U}\_{k,r} < \left(\mathbf{L}\_{k,r+1} - \mathbf{U}\_{k,r}\right) \times \left(1 - b\_{k,r}\right) \forall r \forall k \tag{22}$$

$$I\_k^I = \sum\_{r=1}^K b\_{k,r} \times \mathbb{S}\_{k,r} \forall k \tag{23}$$

Although several objectives are considered simultaneously in this work, it is important to note that these objectives conflict with each other. To address such conflict, this work employs a multi-objective approach known as fuzzy optimization. Fuzzy optimization is a simple multi-objective optimization approach that was founded upon the fuzzy decision-making approach introduced by Bellman and Zadeh [40]. Zimmermann [41] then extended the fuzzy approach to deal with linear and non-linear programming problems that contain multiple objectives. In this approach, a continuous interdependence variable, λ, which is also known as the degree of satisfaction, is introduced. Every fuzzy constraint will be satisfied partially at least to λ. Thus, the multi-objective functions in the optimization can be integrated into a single objective function within the optimization framework. Fuzzy optimization has been widely used to determine the optimum process alternative that considers all aspects simultaneously based on the pre-defined limits. The fuzzy optimization model is expressed as:

$$
\lambda = \sum\_{k=1}^{K} \lambda\_k b\_k \tag{24}
$$

$$\sum\_{k=1}^{K} b\_k = 1\tag{25}$$

where integer variable, *bk* is used to indicate the existence (or absence) of λ*<sup>k</sup>* of process *k*. Note that the formulation in Equation (24) results in the model being non-linear. The optimization objective of λ is maximized, subject to the predefined upper and lower bounds. All flexible targets (GP, HQI, ISI and WAR) are predefined as fuzzy goals which are given by a linear membership function bounded by the upper (GPU, HQIU, ISIU, WARU) and lower limits (GPL, HQIL, ISIL, WARL). λ of each processes *k* is determined using Equations (26)–(29):

$$\frac{\mathrm{GP}\_k - \mathrm{GP}^\mathrm{L}}{\mathrm{GP}^\mathrm{U} - \mathrm{GP}^\mathrm{L}} \ge \lambda\_k \qquad \forall k \tag{26}$$

$$\frac{\text{WAR}^{\text{U}} - \text{WAR}\_{k}}{\text{WAR}^{\text{U}} - \text{WAR}^{\text{L}}} \ge \lambda\_{k} \qquad \forall k \tag{27}$$

$$\lambda \frac{\text{ISI}^{\text{U}} - \text{ISI}\_{k}}{\text{ISI}^{\text{U}} - \text{ISI}^{\text{L}}} \ge \lambda\_{k} \tag{28}$$

$$\frac{\mathrm{HQI}^{\mathrm{U}} - \mathrm{HQI}\_{k}}{\mathrm{HQI}^{\mathrm{U}} - \mathrm{HQI}^{\mathrm{L}}} \geq \lambda\_{k} \qquad \forall k \tag{29}$$

Following this, the model is solved via fuzzy optimization to select a more sustainable process configuration based on the aforementioned four objectives. Based on the optimized solution, the operating capacities, *xj,k* of the selected process configuration are then examined. If the operating capacities do not differ from the values used in the preliminary simulation, then the selected process *k* can be recommended for the next process design phase. However, if the operating capacities do differ from the preliminary simulation, the selected process *k* would be re-simulated in accordance to the optimized operating capacities. Once the selected process *k* is re-simulated, the operating capacities are checked. In the event where the operating capacities in the simulation exceed an allowable error margin (i.e., ≥10%), the preliminary simulation must be revisited for troubleshooting. For cases where the operating capacities are within the allowed error margin, the re-simulated conditions for process *k* can be recommended for the next process design phase. Note that the error margin can be changed based on several factors such as type of industry, experience and knowledge of the decision-maker, and type of equipment used.

#### **4. Case Study**

To illustrate the proposed framework, a case study on 1,4-butanediol (C4H10O2) production process is presented. 1,4-butanediol is a colorless and non-corrosive solution which is widely used as a solvent in the industry to manufacture elastic fibers, technical plastics and polyurethanes [42]. Based on literature review, Reppe [42] and Davy [43,44] processes are the two common and well-established alternative processes available to produce 1,4-butanediol. Table 1 shows the reaction pathway for both processes.

**Table 1.** Reaction pathway for the production of 1,4-Butanediol (C4H10O2).


As shown in Table 1, the raw materials needed for Reppe process are formaldehyde (C2H2), acetylene (CH2O) and hydrogen (H2). Formaldehyde and acetylene are fed into two batch reactors which operate at 5 bar, 80 ◦C and are arranged in parallel. Formaldehyde is reacted with acetylene in the presence of copper, bismuth and silicon dioxide supporting catalyst to produce 1,4-butynediol as an intermediate product with a conversion rate of 60%. Aqueous sodium hydroxide and hydrochloric acid are used to maintain the pH of the medium in the reactor within the range of pH 6 to 8. Based on the reaction, process synthesis is started by establishing the remaining unit operations in the process. For example, the product from the batch reactors is sent to a buffer tank and is stored for approximately 1–30 h before feeding to a distillation column for further separation. After separation, the bottom product stream of the distillation column is rich in 1,4-butynediol and is fed to a trickling and hydrogenation reactor for further reaction to produce 1,4-butanediol as the main product with a conversion rate of 90%. The operating condition of this reactor is 160 ◦C and 150 bar. Figure 3 shows the process flow diagram for the Reppe process.

For the Davy process, the raw materials are dimethyl maleate (DMM), C6H8O4 and hydrogen, H2. DMM is completely vaporised with H2 before feeding into the first fixed bed reactor. The fixed bed reactor operates at 15 bar, 150 ◦C and consists of two bed of catalysts which are palladium on alumina and copper-zinc oxide. In the reactor, DMM is fully converted into dimethyl succinate (DMS) through hydrogenation process in the presence of palladium on alumina as catalyst in the first reaction stage. In the second reaction stage, DMS is mainly converted into gammabutyrolactone (GBL) with trace amounts of tetrahydrofuran (THF) and 1,4-butanediol in the presence of copper-zinc oxide as catalyst. The conversion rate of the second reaction is 97%. Based on these two reactions, process synthesis is conducted to list all the other required unit operations. For instance, the vapour-phase product from the reactor is condensed at 30 ◦C and is fed into a gas-liquid separator to remove the unreacted H2. The unreacted H2 is then being recycled to the fixed bed reactor for further reaction. Then, the GBL rich stream from the first reactor is fed into the second hydrogenation reactor for further reaction to produce 1,4-butanediol with trace amounts of methanol and THF as by-products. The operating condition of the second reactor is 210 ◦C and 75 bar. The catalyst used in the second hydrogenation reactor is copper-zinc oxide and the conversion of GBL to 1,4-butanediol is around 95%. The PFD for the Davy process is shown in Figure 4. The annual production rates for both Reppe and Davy process are assumed to be 60,000 tonnes of 1,4-butanediol, operating at 8000 h/y. The acceptable range of fluctuation on the amount of the product output given by the decision maker is within ±30%. The average height of the main unit operations in the process is assumed to be below 7 m [45]. Wind speed is assumed to be 4 m/s, which is a typical value for outdoor facilities, since local average wind speed is not available [45].

Based on the information mentioned above, process simulation can be performed to analyse the performance for both Reppe and Davy processes. In this work, both processes are simulated with the aid of commercial process simulation software, Aspen HYSYS version 8.8 (Aspentech, Bedford, MA, USA, 2014) [46] (see simulation flowsheets in Figures S1,S2 and simulation files in Supplementary Materials respectively). Settings and parameters used for both simulations are summarized in Tables S1 and S2 (Supplementary Materials) respectively.

Based on the simulation result, data such as input and output mass flow rates and energy consumption are then extracted. The extracted data was then used to develop the input-output models for both processes (see Tables S3 and S4 in Supplementary Materials).

In the following step, the input-output models are evaluated based on economic, environmental, health and safety aspects. For the environmental aspects, the weighting factor of the impact at category *l,*α*<sup>l</sup>* is assumed as one. This indicates that the impact for each category is equally important. Meanwhile, Table 2 shows the information obtained for both Reppe and Davy processes prior to the optimisation step. As shown in Table 2, the Davy process has higher initial GP and WAR; however, with a lower score of HQI and ISI, compared with the Reppe process. The latter has a lower initial GP mainly due to its higher cost of raw materials. Its higher HQI score is mainly due to the presence of harmful chemical components in the process, such as formaldehyde and 1,4-butynediol. Note that the threshold limit (TLV) of formaldehyde and 1,4-butanediol are very low, that is, 1.228 mg/m<sup>3</sup> and 0.5 mg/m3 respectively. Besides, the higher ISI score for the Reppe process is mainly due to its highly exothermic side reaction involving formaldehyde, acetylene and hydrogen to form 2-propanol. Note that the Reppe process has a lower score of WAR (i.e., it is safer) because it converts the hazardous raw material (formaldehyde) to a product (1,4-butanediol) with a lower environmental impact. As such, the selection between these two processes is complex, particularly when four different aspects are considered simultaneously.

**3.**Processflow diagramoftheReppe

**Figure 4.** Process flow diagram of the Davy Process.


**Table 2.** Data obtained for Reppe and Davy processes before optimisation.

In order to determine the optimum process with consideration of all aspects simultaneously, fuzzy optimisation is used. Based on Equations (26)–(29), it is noted that lower and upper limits for GP, HQI, ISI and WAR are required. To obtain these values, *GPk* is first maximized to determine the upper limit GPU. The corresponding values for *HQIk*, *ISIk* and *WARk* are used as upper limits in HQIU, ISIU and WARU respectively. Following this, *WARk* is minimized to determine WARL. The corresponding *GPk*, *ISIk* and *WARk* values were taken as lower limits GPL, HQIL and ISI<sup>L</sup> respectively. The obtained values for GPU, GPL, HQIU, HQIL, ISIU, ISIL, WARU and WAR<sup>L</sup> are 185.112 <sup>×</sup> <sup>10</sup><sup>6</sup> USD/y, 22.817 <sup>×</sup> 106 USD/y, 0.593, 0.017, 36, 27, 50.519 and 1.045 respectively. The model is then solved by maximizing λ in Equation (24) with constraints in Equations (1)–(23) and (26)–(29). The model was solved using LINGO version 14.0 (Lindo Systems, Chicago, IL, USA, 2015) on a Lenovo P700 with 8 GB RAM and Intel®Core™ i7, 2.60 GHz Processor. The mixed integer non-linear programming (MINLP) model consists of 187 variables, 6 integers and 341 constraints.

Based on the results obtained, the Davy process was selected as it has a higher value of λ*,* that is, 0.467, as compared to that of the Reppe process (0.111). The production of 1,4-butanediol after optimisation is scaled down from 7.9 ton/h to 7.64 ton/h, which is within the acceptable range for the amount of product output. Table 3 shows the comparison of the four aspects before and after optimisation for the Davy process.


**Table 3.** Comparison of the four aspects for the Davy process before and after optimisation step.

After optimisation, it is observed that the operating capacity of the process has been scaled down by 3.3% from the target production to obtain a trade-off among economic, health, safety and environmental aspects. The reason for this is that, in order to achieve the trade-off solution shown in Table 3, the operating capacity for the Davy process has to be scaled down by 3.3%. From here, the new scaled down operating capacity of the Davy process is re-simulated to study the deviations. The net output of the process streams obtained from the developed model has a slight difference to the solution generated from process simulator with a deviation of less than 10%. This is mainly due to the different models used in generating the solution. The solution obtained by the mathematical model is linear, while that used by the process simulator is in rigorous mode. Since the deviation obtained is small, the selected Davy process with its optimized operating capacity can be recommended for the next stage of design.

#### **5. Conclusions**

In conclusion, this work presented a systematic framework to screen and select a sustainable chemical process at the preliminary design stage. The presented framework consists of three main tasks. First, process simulation was carried out to analyse the performance of each alternative process. Based on the simulation, process data such as mass flow rates and energy requirements for each process

unit were determined. These data were then used in the next task, which is to develop an input-output model for each process. Each input-output model was formulated to address four conflicting objectives. These four objectives, namely economic performance, health, safety and environment aspects were assessed using gross profit (GP), the Health Quotient Index (HQI), the Inherent Safety Index (ISI) and the WAR algorithm respectively. Following this, fuzzy optimization was used in the third task to optimize the abovementioned aspects simultaneously and to select the most sustainable process. A sustainable process in this context refers to the process with the highest degree of satisfaction (λ) among all four aspects considered. The proposed framework was illustrated with a case study on 1,4-butanediol production, where the proposed framework was used to determine the most sustainable process to produce 1,4-butanediol. Two processes, namely the Reppe and Davy processes were considered. Each process was simulated using Aspen HYSYS Version 8.8 and subsequently underwent input-output modelling. The developed input-output models were then optimized using fuzzy optimization. Results from the case study indicate that the Davy process is more sustainable compared to the Reppe process. In particular, the Davy process had the highest λ value. Furthermore, the results suggest that the operating capacity of the Davy process should be scaled down by 3.3% to meet the trade-off scores determined by the λ value obtained. Based on this recommendation, the Davy process was re-simulated. It was found that there was a small percentage of deviations between the mathematical and process simulation models. This proved that the developed model is feasible for determining a more sustainable process configuration. For future work, a selection of different unit operations may be incorporated into the model to determine their impact on flowsheet selection.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2227-9717/7/4/200/s1, Figure S1: Simulation Flowsheet of Davy Process via Aspen HYSYS version 8.8; Figure S2: Simulation Flowsheet of Reppe Process via Aspen HYSYS version 8.8; Table S1: Simulation Settings for Davy Process; Table S2: Simulation Settings for Reppe Process; Table S3: Input-Output Table for Reppe Process; Table S4: Input-Output Table for Davy Process; Simulation File for Davy Process (Filename: Davy Process); Simulation File for Reppe Process (Filename: Reppe Process).

**Author Contributions:** Conceptualization, D.C.Y.F. and D.K.S.N.; methodology, S.Y.T., K.B.C., D.C.Y.F., D.K.S.N. and M.H.H.; software, S.Y.T., K.B.C., B.H.H., A.J.W.L. and V.A.; validation, V.A. and D.K.S.N.; formal analysis, S.Y.T., K.B.C. and B.H.H.; investigation, S.Y.T., K.B.C., B.H.H.; resources, D.K.S.N.; data curation, B.H.H., A.J.W.L.; writing—original draft preparation, S.Y.T., V.A.; writing—review and editing, V.A., D.C.Y.F., M.H.H. and D.K.S.N.; visualization, S.Y.T.; supervision, V.A., D.C.Y.F., M.H.H. and D.K.S.N.

**Funding:** This research was funded by Ministry of Higher Education, Malaysia through the LRGS Grant (LRGS/2013/UKM-UNMC/PT/05).

**Acknowledgments:** The authors would like to acknowledge the financial support provided from the Ministry of Higher Education, Malaysia through the LRGS Grant (LRGS/2013/UKM-UNMC/PT/05).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Nomenclature**





#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Symmetry Detection for Quadratic Optimization Using Binary Layered Graphs**

#### **Georgia Kouyialis †, Xiaoyu Wang and Ruth Misener \***

Department of Computing, Imperial College London, London SW7 2AZ, UK;

g.kouyialis14@imperial.ac.uk (G.K.); xiao.wang10@imperial.ac.uk (X.W.)

**\*** Correspondence: r.misener@imperial.ac.uk; Tel.: +44-207-594-8315

† Current address: Schlumberger Cambridge Research.

Received: 16 May 2019; Accepted: 5 November 2019; Published: 9 November 2019

**Abstract:** Symmetry in mathematical optimization may create multiple, equivalent solutions. In nonconvex optimization, symmetry can negatively affect algorithm performance, e.g., of branch-and-bound when symmetry induces many equivalent branches. This paper develops detection methods for symmetry groups in quadratically-constrained quadratic optimization problems. Representing the optimization problem with adjacency matrices, we use graph theory to transform the adjacency matrices into binary layered graphs. We enter the binary layered graphs into the software package nauty that generates important symmetric properties of the original problem. Symmetry pattern knowledge motivates a discretization pattern that we use to reduce computation time for an approximation of the point packing problem. This paper highlights the importance of detecting and classifying symmetry and shows that knowledge of this symmetry enables quick approximation of a highly symmetric optimization problem.

**Keywords:** symmetry; quadratic optimization; quadratically-constrained quadratic optimization

#### **1. Introduction**

When the optimization variables can be permuted without changing the structure of the underlying optimization problem, we say that the *formulation group* of an optimization problem is *symmetric* [1,2]. For motivation, consider the circle packing problem illustrated in Figure 1 [3]. Given an integer *n* > 0, the circle packing problem asks: what is the largest radius *r* for which *n* non-overlapping circles can be placed in the unit square? Costa et al. [3] show that the formulation group, i.e., a subgroup of symmetry group generated by permuting variables and constraints, is isomorphic to a symmetry group created by permuting the variable indices and switching the two coordinates in a unit square (*C*<sup>2</sup> × *Sn*). Solution methods for nonconvex optimization problems lacking symmetry-aware formulations and/or solution procedures may end up exploring all of these equivalent solutions. In other words, symmetry may cause classical optimization methods such as branch-and-bound to explore many unnecessary subtrees.

More generally, a number of authors have considered a range of symmetry detection methods, e.g., for constraint programming [4], integer programming [1,5–8], and mixed-integer nonlinear optimization [2]. These automatic symmetry detection methods can then be used to mitigate the computational difficulties caused by symmetries, e.g., with symmetry-breaking constraints [9–11], objective perturbation [12], specialized branching strategies [13,14], cutting planes [15,16], and extended formulations [17]. The recent computational comparison of Pfetsch and Rehn [18] indicates that these state-of-the-art symmetry handling methods expedite the solution process for the MIPLIB 2010 instances and additionally enable more instances to be solved in a time limit.

**Figure 1.** Given an integer *n* > 0, the circle packing problem asks: what is the largest radius *r* for which *n* non-overlapping circles can be placed in the unit square? Already for *n* = 2, there are four equivalent solutions [3]; these solutions are related to one another via rotations and reflections.

Researchers have also developed symmetry-handling methods for specific applications including covering design [19], circle packing [3,20], scheduling [21], transmission switching [22], unit commitment [23–26], and heat exchanger network synthesis [27,28]. As a concrete example of the type of contributions researchers have made, consider a job shop scheduling problem that minimizes makespan on two identical machines. Good scheduling formulations and/or solution procedures, e.g., Maravelias and Grossmann [29], Maravelias [30], and Mistry et al. [31], will implicitly exclude two of the three equivalent solutions diagrammed in Figure 2.

**Figure 2.** To observe symmetries that may arise in scheduling, consider a job shop scheduling problem that minimizes makespan on two identical machines. Lacking symmetry-aware formulations and/or solution procedures, a solution procedure may end up exploring all three of these equivalent solutions.

This paper develops detection methods for symmetry groups in quadratically-constrained quadratic optimization problems. Representing the optimization problem with adjacency matrices, we use graph theory to transform the adjacency matrices into binary layered graphs. We enter the binary layered graphs into the software package nauty that generates important symmetric properties of the original problem. Symmetry pattern knowledge motivates a discretization pattern that we use to reduce computation time for an approximation of the point packing problem.

#### **2. Formulation Symmetry for Quadratically-Constrained Quadratic Optimization Problems**

Consider the quadratically-constrained quadratic optimization problem (QCQP):

$$\begin{aligned} \min\_{\mathbf{x} \in \mathbb{R}^n} \quad & f\_0(\mathbf{x}) \\ \text{s.t.} \quad & f\_k(\mathbf{x}) \le 0 \qquad \forall k = 1, \dots, m \\ & \mathbf{x}\_i \in \left[ \mathbf{x}\_i^L, \mathbf{x}\_i^{\mathrm{ul}} \right] \quad \forall i = 1, \dots, n \end{aligned} \tag{\text{QCQP}}$$

where:

$$f\_k(\mathbf{x}) = \sum\_{i=1}^n \sum\_{j=1}^n a\_{ij}^k \mathbf{x}\_i \mathbf{x}\_j + \sum\_{i=1}^n a\_{i0}^k \mathbf{x}\_i + a\_{00}^k \forall k = 0, \dots, m,\tag{1}$$

with finite variable bounds *x<sup>L</sup> <sup>i</sup>* , *<sup>x</sup><sup>U</sup> <sup>i</sup>* <sup>∈</sup> <sup>R</sup>, <sup>∀</sup>*<sup>i</sup>* and coefficients *<sup>α</sup><sup>k</sup> ij* <sup>∈</sup> <sup>R</sup> for *<sup>i</sup>*, *<sup>j</sup>* ∈ {0, . . . , *<sup>n</sup>*}, *<sup>k</sup>* ∈ {0, . . . , *<sup>m</sup>*}.

To represent symmetry in the QCQP formulation, consider *Sn*, the symmetric group of order *n* formed by the *n*! possible permutation operations. The *formulation group* of QCQP , or G*QCQP*, is the set of variable index permutations that preserve the objective and constraint structure [2]. For a variable index permutation *π* ∈ *Sn*, we seek the constraint index permutations *σ* ∈ *Sn* that maintain both the objective value and the constraint structure on the feasible domain *dom*(*f*) where *f* = [ *f*1, ... , *fm*]. More formally:

**Definition 1** (Formulation group of QCQP)**.**

$$\mathcal{G}\_{\mathsf{QCQP}} = \left\{ \pi \in \mathcal{S}\_{\mathsf{n}} \mid \forall \mathbf{x} \in dom(f\_0) \, f\_0(\pi \mathbf{x}) = f\_0(\mathbf{x}) \land \forall \mathbf{x} \in dom(f) \, \exists \sigma \in \mathcal{S}\_{\mathsf{m}} \left( \sigma f(\pi \mathbf{x}) = f(\mathbf{x}) \right) \right\}.$$

Because *dom*(*f*) may be nonconvex and difficult to compute, this paper considers a <sup>G</sup>*QCQP* restriction that enforces symmetry on the entire box bounds, i.e., we assume that *dom*(*f*) = *xL*, *xU* for the purpose of computing G*QCQP*. The next subsections represent formulation symmetry in two ways: (i) expression graphs in Section 2.1 and (ii) tensors in Section 2.2. Representing formulation symmetry using expression graphs is due to Liberti [2] and the tensor representation is new to this paper.

#### *2.1. Symmetry Detection with Expression Graphs*

One option to compare two functions is to compare their expression trees, i.e., a directed acyclic graph representation of each function that incorporates the relevant operations, constants, and variables [2]. These expression tree models were first developed for mixed integer nonlinear optimization (MINLP) by Smith and Pantelides [32] and are common in most global MINLP solvers [33–39] and other MINLP-related software [40–42]. Figure 3 illustrates a simple example of an expression tree for 3*x*<sup>1</sup> + 2*x*<sup>2</sup> <sup>4</sup> + 2*x*2*x*3. A tree comparison algorithm may recursively compare two trees to determine equivalence [2]. More advanced implementations may detect equivalent but differently-formulated expressions, e.g., (*x*<sup>1</sup> + *x*2)<sup>2</sup> versus *x*<sup>2</sup> <sup>1</sup> + <sup>2</sup>*x*1*x*<sup>2</sup> + *<sup>x</sup>*<sup>2</sup> 2.

**Figure 3.** Example of an expression tree for 3*x*<sup>1</sup> + 2*x*<sup>2</sup> <sup>4</sup> + 2*x*2*x*3.

With a directed acyclic graph representation, Liberti [2] computes the formulation symmetry group using the graph isomorphism problem, i.e., a problem that can be solved using off-the-shelf software nauty [43]. Liberti [2] also proves how to map the automorphism group of a directed acyclic graph to the formulation group of the original MINLP.

#### *2.2. Symmetry Detection with Tensors*

As an alternative to the expression tree representation, Figure 4 illustrates that QCQP can be represented as a tensor: *<sup>A</sup>QCQP* <sup>∈</sup> <sup>R</sup>(*n*+1)×(*n*+1)×(*m*+2*n*). Each of the two dimensions (*<sup>n</sup>* <sup>+</sup> <sup>1</sup>) corresponds to a constant term and the variables. Each two-dimensional slice of the tensor corresponds to the constant, linear, and quadratic terms in a constraint. The first *m* slices correspond to Equation (1) and have entries *a<sup>k</sup> ij*. The next 2*<sup>n</sup>* slices correspond to the box constraints, i.e., *xi* ≥ *<sup>x</sup><sup>L</sup> <sup>i</sup>* and *xi* ≤ *<sup>x</sup><sup>U</sup> <sup>i</sup>* , ∀*i*.

The formulation group of this representation is:

$$\mathcal{G}\_{\mathrm{QCD},T} = \left\{ \pi \in \mathcal{S}\_{\mathbb{R}} \mid \forall \mathbf{x} \in dom(f\mathfrak{d}) \left( f\mathfrak{d}(\pi \mathbf{x}) = f\mathfrak{d}(\mathbf{x}) \land \forall \mathbf{x} \in dom(f) \; \exists \sigma \in \mathcal{S}\_{\mathbb{R}} \left( \mathcal{A}\_{\mathrm{QCD}}(\pi, \pi, \sigma) = \mathcal{A}\_{\mathrm{QCD}} \right) \right) \right\}.$$

**Figure 4.** Tensor representation of the symmetry.

#### 2.2.1. Sparse Tensor Representation

For a given tensor *AQCQP*, consider a sparse representation, illustrated in Figure 5, that reduces the memory required to store the tensor. Instead of storing the entire tensor, we store arrays of length *s* where *s* is the number of nonzero entries in QCQP. The first array, *M* <sup>=</sup> (*M*1, ... , *Ms*) stores all nonzero entries *α<sup>k</sup> ij* of QCQP. The next three arrays, *<sup>I</sup>* = (*I*1, ... , *Is*), *<sup>J</sup>* = (*J*1, ... , *Js*), *<sup>K</sup>* = (*K*1, ... , *Ks*) represent the indices corresponding to the nonzero *α<sup>k</sup> ij* entries. The maximum size of *<sup>s</sup>* is (*<sup>n</sup>* + <sup>1</sup>)2(*<sup>m</sup>* + <sup>2</sup>*n*), but, in practice, most arrays will be significantly shorter.

**Figure 5.** Sparse tensor representation models *AQCQP* as 4 arrays with the nonzero entry *α<sup>k</sup> ij* in *<sup>M</sup>* and arrays *I*, *J*, *K* holding the index.

#### 2.2.2. Converting Matrices to Edge-Labeled, Vertex-Colored Graphs

We convert the sparse tensor representation of *AQCQP* into an edge-labeled, vertex-colored graph. Given the edge-labeled, vertex-colored graph, generating graph automorphisms to the original problem symmetries is well-known [1,5,44,45]. To construct the edge-labeled, vertex-colored graph, consider a graph *<sup>G</sup>* = (*V*, *<sup>E</sup>*, *<sup>c</sup>*) corresponding to an instance *<sup>M</sup>*, *<sup>I</sup>*, *<sup>J</sup>*, *<sup>K</sup>*. The function *<sup>c</sup>* : *<sup>E</sup>* <sup>→</sup> *<sup>r</sup>*, for *r* ∈ {0, ... , - − 1} is an edge coloring where - <sup>∈</sup> <sup>Z</sup><sup>+</sup> is the number of different coefficients in *<sup>M</sup>*. Each unique element in *<sup>M</sup>* is stored in a vector *<sup>U</sup>* <sup>∈</sup> *<sup>R</sup>*-. We also partition (color) the vertex set into four subsets: a set *VF* representing the objective function, *VC* nodes for the constraints, a constant node *VS*, and *VR* variable nodes. The automorphism definition prevents vertices from being mapped onto a vertex of a different color, so these colors prevent, for example variables becoming constraints. The equivalence relation is [2]:

$$\begin{aligned} \forall \iota, v \in V \flat \iota \sim v &\implies \left(\iota, v \in V\_{\mathbb{F}} \wedge \ell(\iota) = \ell(v)\right) \vee \left(\iota, v \in V\_{\mathbb{C}} \wedge \ell(\iota) = \ell(v)\right) \\ &\lor \left(\iota, v \in V\_{\mathbb{S}} \wedge \ell(u) = \ell(v)\right) \vee \left(\iota, v \in V\_{\mathbb{R}} \wedge \ell(u) = \ell(v)\right) .\end{aligned}$$

Figure 6 illustrates the edge-labeled, vertex-colored graph. Initially, the edge set is empty *E* = ∅. For *i* = {0, ... ,*s*} where *s* = |*M*|, add an edge *v* (*r*) *Ii* to *v* (*r*) *Ki* , i.e., from a vertex in the set that represents the constant element / variables to a vertex in the set of the objective function / constraints, with the relevant color. The graph construction incorporates edges between variable nodes *VR* for the quadratic bilinear terms. For *i* = {0, . . . ,*s*}:


**Figure 6.** The tensor *AQCQP* as an edge-labeled, vertex-colored graph.

#### **3. Formulation Symmetry Detection via Binary Layered Graphs**

The software nauty [43], which detects symmetry, accepts vertex-colored graphs but does not accept the Section 2.2.2 edge-labeled, vertex-colored graphs. Thus, we associate edge colors with layers in a graph and transform the edge-labeled, vertex-colored graph into a vertex-colored graph. Since the transformation from an edge-labeled, vertex-colored graph to a vertex-colored graph is isomorphic [43], the transformation does not lose anything. Using the resulting *binary layered graphs*, we generate the automorphism group and find symmetry in the original QCQP.

To convert a graph *G* = (*V*, *E*, *c*) with colors into an - *- layered graph* [43], we replace each vertex *vj* ∈ *V* with a fixed connected graph of vertices *v* (0) *<sup>j</sup>* , ... , *v* (-−1) *<sup>j</sup>* . If an edge (*vj*, *vj*) has color *r*, add an edge from *v* (*r*) *<sup>j</sup>* to *v* (*r*) *<sup>j</sup>* . Finally, we partition the vertices by the superscripts, *Vr* = {*v* (*r*) <sup>0</sup> , ... , *v* (*r*) *<sup>n</sup>*−1}. Alernatively, a binary representation avoids too many layers in *G* when the number of colors is large.

**Definition 2** (Binary Layered Graph)**.** *Let* - <sup>∈</sup> <sup>Z</sup><sup>+</sup> *be the number of edge labels of G. A binary layered graph is a vertex-colored graph where the number of layers L* = log2 (-+ 1) *matches a binary representation.*

We assign a unique positive integer *<sup>μ</sup>*(*z*) to each *<sup>z</sup>* <sup>∈</sup> *<sup>U</sup>* and map edge labels *<sup>μ</sup>*(*z*) to a binary representation that switches on/off parameters *ct* to represent the edge colors as layers. If *ct* = 1, add a new edge from *v<sup>t</sup> <sup>i</sup>* to *<sup>v</sup><sup>t</sup> <sup>j</sup>* for every *ct* ∈ {*c*1,..., *cL*−1}:

$$\mu(z) = 2^{L-1} \cdot c\_{L-1}(z) + 2^{L-2} \cdot c\_{L-2}(z) + \dots + 2^0 \cdot c\_0(z), \text{ for } c\_t \in \{0, 1\}, t = \{0, \dots, L-1\}. \tag{2}$$

Figure 7 illustrates the resulting binary labeled graph with its *L* = log2 (- + 1) + 2 layers. There are vertices for the objective function and each constraint and layers of copies of these constraints (connected with vertical edges). The horizontal edges encode the problem coefficients. On the upper part of Figure 7, there are vertices for a constant element and each variable and a layer of variable copies (connected with vertical edges). Here, the horizontal edges and loops distinguish the linear and bilinear terms. Algorithms 1 and 2 summarize computing the vertex and edge sets, respectively.

**Figure 7.** Binary layered graph representation of QCQP using the tensor representation *AQCQP*.

**Algorithm 1** Algorithm constructing the vertex set


**Algorithm 2** Algorithm constructing the edge set

```
1: procedure G=(V, E)
2: V ← V
3: E ← ∅
4: for s = 0 → L − 4 do  Vertical edges between copies of vertices
5: for k = 0... m do  Copies of vertices representing the constraints
6: E = E ∪ E ∪ (v
                  (s)
                  k , v
                     (s+1)
                     k )
7: end for
8: end for
9:  Vertical edges between copies of vertices
10: for j = 0... n do  Copies of vertices representing the variables
11: E = E ∪ (v
             (L−2)
             j , v
                  (L−1)
                  j )
12: end for
13: for F = 0, . . . N − 1 do  Bilinear terms
14: if IM(F) = JM(F) then  Add a loop
15: E = E ∪ (v
               (L−1)
                IM(F)
                   , v
                    (L−1)
                    IM(F)
                        )
16: else
17: if IM(F) < JM(F) then  Add an edge
18: E = E ∪ (v
                  (L−1)
                  IM(F)
                     , v
                      (L−1)
                       JM(F)
                          )
19: else
20: E = E ∪ ∅
21: end if
22: end if
23: end for
24: for F = 0... N − 1 do
25: for k = 0... m do
26: if KM(F) = k then
27: if IM(F) = 0 then
28: E = E ∪ (v
                    (s)
                    k , v
                       (L−2)
                       JM(F)
                          )
29: else
30: E = E ∪ (v
                    (L−1)
                    IM(F)
                       , v
                         (L−1)
                         JM(F)
                            )
31: end if
32: end if
33: end for
34: end for
35: return G
36: end procedure
```
#### **4. Numerical Discussion and Comparison to the State-of-the-Art**

The following example incorporates the algorithms proposed in this paper. We construct the binary labeled graph and then enter it into nauty through the dreadnaut command line interface:

$$\max\_{\mathbf{x}\_1, \mathbf{x}\_2, \mathbf{x}\_3, \mathbf{x}\_4 \in [0, 1]} \quad 3\mathbf{x}\_1 + 3\mathbf{x}\_4 + 2\mathbf{x}\_2 \mathbf{x}\_3 \tag{c\_0},$$

$$\begin{aligned} \mathbf{x}\_2 + \mathbf{x}\_1^2 + 1 &\le 0 \\ \mathbf{x}\_3 + \mathbf{x}\_4^2 + 1 &\le 0 \\ \mathbf{x}\_2 + \mathbf{x}\_3 + 1 &\le 0 \end{aligned} \qquad \begin{aligned} (c\_0), \\ (c\_1), \\ (c\_2), \\ (c\_3). \end{aligned}$$

The optimization problem has sparse matrix representation: *M* = (3, 3, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1), *I* = (0, 0, 2, 0, 0, 1, 0, 0, 4, 0, 0, 0), *J* = (1, 4, 3, 0, 2, 1, 0, 3, 4, 0, 2, 3), *K* = (0, 0, 0, 1, 1, 1, 2, 2, 2, 3, 3, 3), vector of unique elements *<sup>U</sup>* = [1, 2, 3], and *<sup>L</sup>* <sup>=</sup> log2 <sup>4</sup> <sup>=</sup> 2 layers.

Equation (2) computes the binary representation of each unique element, e.g., 3 = 2<sup>1</sup> + 20 indicates that there is an edge between vertices on layer zero and another edge between the same vertices on layer 1. The graph consists of four layers and |*V*| = 18, one associated with a constant element and one with the objective function and the rest for the problem variables and constraints. The left-hand side of Figure 8 illustrates the graph representation. Nauty generates permutations: *π* = (1, 2)(5, 6)(9, 12)(10, 11)(14, 17)(15, 16). To see how these Nauty-generated permutations usefully explain the symmetry properties of the entire problem, observe: (i) Permutations (1, 2)(5, 6), as shown in Figure 8, permute the constraints *c*1, *c*<sup>2</sup> and (ii) Permutations (9, 12)(10, 11) are associated with the variables *x*1, *x*<sup>4</sup> and *x*2, *x*<sup>3</sup> with (14, 17)(15, 16) their copies. These permutations therefore allow us to automatically calculate the formulation group G = (*x*1*x*4)(*x*2*x*3).

**Figure 8.** Illustration of the example problem using the binary labeled graph representation (**left**) and the directed acyclic graph representation (**right**).

The right-hand side of Figure 8 uses Section 2.1 to develop a directed acyclic graph representation for the same problem. The graph colors represent the vertex partitioning that enables node exchanges. In this case, the directed acyclic graph representation uses a smaller number of vertices and edges than the tensor-based representation. However, the representations generate the same formulation group.

**Comparison.** To evaluate the trade-offs between the tensor and directed acyclic graph representations, observe that both methods will search for the same formulation group symmetries. However, the tensor representation may be especially useful when working with problems with many differently-valued coefficients, i.e., the logarithmic number of layers may reduce the number of nodes. The function assigning integer values to the problem coefficients lets us work not only with 0–1 coefficients, but also with any other value.

#### **5. Exploiting Symmetry in the Point Packing Problem**

Once symmetry has been detected, we can use our knowledge of the symmetry to mitigate the computational difficulties caused by symmetry. Here, we focus on solving the point packing problem [46–49]. The point packing problem concerns packing *n* points to a unit square. The aim is to maximize the in-between distance between any two points:

$$\begin{aligned} \max \quad & \quad \theta\\ \text{subject to} \quad & (x\_i - x\_j)^2 + (y\_i - y\_j)^2 \ge \theta & 1 \le i < j \le n\_i\\ & 0 \le x\_i \le 1, 0 \le y\_i \le 1 & 1 \le i \le n\_i \end{aligned}$$

where *θ* denotes the minimum distance between any two points. To approximate this problem, consider a grid approach that approximates the optimal solution by adding grid lines to the unit square and forcing the points to be placed only to the vertices generated by these grid lines. Note that the point packing problem has significant applications, e.g., in placing mobile phone towers.

For *n* points, there will be at most a *k* ∗ ×*k*∗ grid, where *k*∗ is the smallest number whose square is the least integer that is greater than *<sup>n</sup>*, i.e., (*<sup>k</sup>* ∗ −1)<sup>2</sup> < *<sup>n</sup>* and *<sup>k</sup>*∗<sup>2</sup> ≥ *<sup>n</sup>*. In other words, we add at most 2*k*∗ grid lines. On this *k* ∗ ×*k*∗ grid, points will occupy most of the vertexes and the unit length of spacing has been maximized. This approach, unfortunately, has the potential of missing the optimal solution. Consider fitting six points to the square. The most fitting grid is 3 × 3 and the optimal *θ* we achieve is 0.5. However, the optimal solution of 0.6009 is achieved by the arrangement shown in Figure 9, which is not available on a 3 × 3 grid. Although an approximation, *k*∗ is still a useful pruning tool, e.g., points need to be at least <sup>1</sup> *<sup>k</sup>*∗−<sup>1</sup> away from any other point.

**Figure 9.** Optimal arrangement of six points.

#### *Exhaustive Search and 2D Symmetry Removal*

First, consider Algorithm 3, an exhaustive search method. To break the symmetries, start by calculating how many grids can possibly have points and the upper limit on the number of points on any occupied grid lines. Once calculated, the complete set of all possible combinations on the *x*-axis is calculated. For example, on a 5 × 5 grid, one possible *x*-coordinate setup is (2, 0, 1, 0, 2). If we label the five horizontal grid lines from 1 to 5, this setup means that there are two points on grid 1, no points on grid 2, one point on grid 3, no points on grid 4 and two points on grid 5. To remove *x*-axis symmetries, we consider reflections as a duplication, i.e., only one of (1, 0, 2, 0, 2) and (2, 0, 2, 0, 1) are considered.


**Input:** number of points *n*, number of grids *k*, number of occupied grids *m*. **Output:** The optimal solution *d*∗


Algorithm 3 considers the unique *x*-axis combinations. The *y*-axis, i.e., the horizontal grid lines, should preserve the same characteristics. Thus, in the improved Algorithm 4, we generate the set *O* of possible point arrangement such as (2, 0, 1, 0, 2), and this is applied to both horizontal and vertical grid lines; in other words, we now consider the product *O* × *O* to give the exact coordinates of all *n* points.


**Input:** number of points *n*, number of grids *k*, number of occupied grids *m*. **Output:** The optimal solution *d*∗


**Example 1.** *Consider for five points on a* 5 × 5 *grid where we have O as* {*O*<sup>1</sup> = (1, 0, 2, 0, 2),*O*<sup>2</sup> = (1, 2, 0, 0, 2), *O*<sup>3</sup> = (1, 0, 0, 2, 2),*O*<sup>4</sup> = (2, 1, 0, 0, 2),*O*<sup>5</sup> = (2, 0, 1, 0, 2)}*. The finalized set of full coordinates would contain* 25 *elements where the product of O with itself is taken. One particular point setup generated by this approach is for example, O*<sup>3</sup> × *O*4*. If we call the vertical grids as V*<sup>1</sup> *to V*5*, respectively, and the horizontal grids as H*<sup>1</sup> *to H*5*, respectively, O*<sup>3</sup> × *O*<sup>4</sup> *would mean that we have the following point locations:*


*For O*<sup>3</sup> × *O*4*, we would have five complete point setups, as shown in Figure 10. In the last setup, we have also added the labels for grids to match what we defined earlier. This diagram contains all possible arrangements of points under this particular orbit partitioning.*

**Figure 10.** Five point setups for *O*<sup>3</sup> × *O*4. (**a**) point setup1; (**b**) point setup2; (**c**) point setup3; (**d**) point setup4; (**e**) point setup5.

To improve further, we can implement pruning mechanisms such as using the bound <sup>1</sup> *<sup>k</sup>*∗−<sup>1</sup> to prune the non-optimal setups. In this particular example, we would have pruned all five setups as none of them satisfy the 0.5 bound from the most fitting grid. This means that *O*<sup>3</sup> × *O*<sup>4</sup> is not the optimal point setup for five points on a 5 × 5 grid.

#### **6. Results and Comparisons**

Both algorithms have been implemented and run on the same devices (HP EliteDesk 800 G2 TWR Intel Core i7-6700 3.4 GHz) to provide effective comparisons. Figure 11 shows the experimental results of 2D Symmetry Removal for packing 6 points. We eliminated the result from Exhaustive Search as it is clear that 2D Symmetry Removal outperforms Exhaustive Search in terms of run-time. We notice that *m*, the number of occupied grids (in the diagram, this corresponds to the occupied vertical grids), has an impact on the run-time. To a reasonable extent, the larger the *m*, the more choices we have regarding where we place the points so it in general takes longer time to compute. Although our strategies effectively convert the point packing QCQP into a mixed-integer linear optimization problems, we could have alternatively designed a branch-and-bound algorithm that is symmetry aware.

**Figure 11.** Run-time for six points on different grids with different number of occupied grids. The line *m* = 4 is above the line *m* = 3.

#### **7. Conclusions**

This paper has explored alternative representations for finding symmetry in formulation groups of a quadratically-constrained optimization problem. We also show that, after knowing the symmetry, we can design significantly better methods to solve the optimization problems. The contributions in this paper are relevant to industrial problems that contain a point packing element [50–52].

**Author Contributions:** Conceptualization, G.K. and R.M.; methodology, G.K.; validation, X.W.; writing, G.K., X.W., and R.M.; supervision, R.M.

**Funding:** This work was funded by an Engineering & Physical Sciences Research Council (ESPRC) Research Fellowship to R.M. (Grant No. EP/P016871/1) and an EPSRC DTP to G.K.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **Abbreviations**

The following abbreviations are used in this manuscript:

MINLP Mixed-integer nonlinear optimization

QCQP Quadratically-constrained quadratic optimization

#### **References**


c 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Modeling and Analysis of Coal-Based Lurgi Gasification for LNG and Methanol Coproduction Process**

#### **Jingfang Gu 1, Siyu Yang 1,\* and Antonis Kokossis 2,\***


Received: 23 May 2019; Accepted: 30 September 2019; Published: 2 October 2019

**Abstract:** A coal-based coproduction process of liquefied natural gas (LNG) and methanol (CTLNG-M) is developed and key units are simulated in this paper. The goal is to find improvements of the low-earning coal to synthesis natural gas (CTSNG) process using the same raw material but producing a low-margin, single synthesis natural gas (SNG) product. In the CTLNG-M process, there are two innovative aspects. Firstly, the process can co-generate high value-added products of LNG and methanol, in which CH4 is separated from the syngas to obtain liquefied natural gas (LNG) through a cryogenic separation unit, while the remaining lean-methane syngas is then used for methanol synthesis. Secondly, CO2 separated from the acid gas removal unit is partially reused for methanol synthesis reaction, which consequently increases the carbon element utilization efficiency and reduces the CO2 emission. In this paper, the process is designed with the output products of 642,000 tons/a LNG and 1,367,800 tons/a methanol. The simulation results show that the CTLNG-M process can obtain a carbon utilization efficiency of 39.6%, bringing about a reduction of CO2 emission by 130,000 tons/a compared to the CTSNG process. However, the energy consumption of the new process is increased by 9.3% after detailed analysis of energy consumption. The results indicate that although electricity consumption is higher than that of the conventional CTSNG process, the new CTLNG-M process is still economically feasible. In terms of the economic benefits, the investment is remarkably decreased by 17.8% and an increase in internal rate of return (IRR) by 6% is also achieved, contrasting to the standalone CTSNG process. It is; therefore, considered as a feasible scheme for the efficient utilization of coal by Lurgi gasification technology and production planning for existing CTSNG plants.

**Keywords:** coproduction; Lurgi syngas; cryogenic separation; methanol synthesis; LNG

#### **1. Introduction**

China is the main source of global energy growth as well as the largest energy consumer in the past 20 years [1]. In 2016, China's natural gas production was 148.7 billion Nm3, with a yearly rate increase of 8.5%. Meanwhile, the total imports are 92 billion m3, with an annual growth rate of 27.6%. However, the total yearly gas consumption is 237.3 billion m3, which is 15.3% higher than that of 2015 [2]. If keeping with the same growth rates, the natural gas will be in insufficient supply in the near future. To alleviate such an energy shortage, the Chinese government encourages the build and operate (B&O) development of coal to conduct synthetic natural gas (CTSNG) projects. Therefore, many CTSNG projects have been launched recently and are being run successfully all over the country, as shown in Table 1.


#### **Table 1.** B&O CTSNG projects in China (2018).

However, CTSNG projects face some challenges. Firstly, the market price of synthetic natural gas (SNG) products is not based on its cost structure, nor according to the guidance from a market mechanism. Departments of China offer a rather lower price to the public, and priority is given to civil use, transportation field, etc. Thus, lower economic returns are common to all CTSNG projects. The price of natural gas has fallen sharply since November 2015 in the country after the National Development and Reform Commission issued a report about price adjustment [2]. In Figure 1, it shows that, during 2015 to 2018, the price of SNG is reduced from 2.75 to 1.82 CNY/Nm3. After 2018, the price further decreased to 1.78 CNY/Nm3. In such case, the market price of SNG is 0.97 CNY/ Nm3 lower than that of 2015. Taking Keqi Coal-Based Gas Project Phase I for an example, in 2017, it produced 1.03 billion Nm3 SNG products, yet with a big deficit of 650 million CNY. It can be seen that if natural gas price keeps fluctuating at a low-price level, these CTSNG projects are likely to face severe losses.

**Figure 1.** China's synthesis natural gas (SNG) market price recording (2014–2018).

Secondly, the SNG product is mainly supplied for civil use like urban heating in the winter. Thus, there is a peak–valley difference of natural gas demand between winter and summer. In most of the northern cities of China, the top demand of natural gas in the heating period is from November to March of the following year. In contrast, the demand in the non-heating period remains at a lower level from April to October. According to statistics from Ji [3], the consumption in the heating period is up to 10 times of that in the non-heating period. Taking the winter of 2018 as an example, the gap between supply and demand of natural gas is about 24 billion m3. Since natural gas cannot be stored for a long time, coal-based gas projects are facing production cuts during the non-heating period, which brings huge economic losses.

Further, in many CTSNG projects, a Lurgi gasification technology has been employed to produce qualitative coal-based syngas for synthesis reaction. Major units in the process can be seen in Figure 2, including coal gasification, water–gas shift, acid gas removal, and methanation synthesis units [4,5].

**Figure 2.** Flowchart of the CTSNG process.

However, the hydrogen:carbon ratio of crude syngas from Lurgi gasifier is about 2.7 [6]. According to the requirement of synthesis gas reaction, it is necessary to use water–gas shift technology to increase that ratio to 3.1 for methanation. However, as the Equation (1) of water–gas shift reaction shows, CO2 emission is inevitably increased in that process [7,8].

$$\text{CO} + \text{H}\_2\text{O} (\text{vapor}) = \text{CO}\_2 + \text{H}\_2 \cdot \Delta \text{H} = 41.19 \text{ kJ/mol} \tag{1}$$

The coproduction process alternative is a practical way to address this challenge. As known, traditional CTSNG processes have a single gaseous coal-based natural gas product. However, it is possible that the same raw materials can be converted to various designed products under the coproduction process structure. Till now, there are different studies devoted to comprehensive processing of coal syngas. These studies prove that coproduction systems can improve the resource utilization and energy efficiencies. Some works are also being proven by demonstration projects.

Yi et al. (2017) studied the modelling and optimization theories of coproduction systems. In their studies, the coproduction process can be very flexible for its integration of technologies and raw material distribution. Besides, it is pointed out that systematic design can improve the process performance like better carbon conversion ratio and improved energy saving size [9–11].

Hao et al. (2015) proposed a coproduction process of methanol and electricity with coal and coke oven gas as raw materials. The new system is compared with the process based on CH4/CO2 dry reforming technology, in terms of exergy efficiency, exergy cost, and CO2 emissions. Through the new system, the exergy efficiency can be increased by 7.8%. Besides, the exergy cost can be reduced by 0.88 USD/GJ and the CO2 emission can be reduced by 0.023 kg/MJ [12–14]. Han et al. (2010) introduced a methanol production and integrated gasification combined cycle power generation system using coal and natural gas as fuel. The syngas derived from natural gas and coal is firstly used for methanol synthesis. The unreacted syngas is used in the power plant as fuel. Comparing with the single production system, the coproduction system can save about 10% of fossil fuels [15]. Tu et al. (2015) found that a methanol and electricity coproduction system can obtain the best benefit when the recycle ratio of unreacted gas is assigned with the value between 2.17 and 4.44, with relative energy saving rate and unit energy production approaching an optimum [16]. Huang et al. (2018) introduced a low-energy CO2 capture process after the water–gas shift unit in a poly-generation process. A part of the unreacted syngas is used to generate power. Energy consumption for CO2 capture is 0.7 GJ/t-CO2, bringing a 40.6% reduction compared to that of the coal-to-methanol process [17].

In addition, Bai et al. (2015) studied a poly-generation system of generating methanol and power with the solar thermal gasification of the biomass. The syngas from the biomass gasification is used to produce the methanol via a synthesis reactor. The un-reacted gas is used for the power generation via a combined cycle power unit. The thermodynamic and economic performances of the system are investigated. A portion of the concentrated solar thermal energy can be chemically stored into the syngas. The highest energy efficiency of the poly-generation system is approximately 56.09%, which can achieve the stable utilization of the solar energy and the mitigation of CO2 emission [18].

Many researchers from outside China are also interested in this field. You et al. (2011) studied the optimal distribution of raw materials in different production routes to maximize the benefit of the coproduction process. A superstructure optimization model is formulated as a mixed-integer nonlinear program to determine the optimal process design, and the proposed framework is applied to a comprehensive superstructure of an integrated shale gas for chemical processing, which involves steam cracking of ethane, propane, *n*-butane, and *i*-butane [19,20].

The above studies are mainly based on thermodynamics to reach a higher energy utilization, achieve a reduction on energy consumption, and realize the optimization of reaction conditions, like gas recycle ratio, operating temperature and pressure, etc. However, studies are less focused on matching products proposal and syngas component ratio, like (H2 − CO2)/(CO + CO2) ratio, which is specified for chemical synthesis.

Considering all difficulties that existing B&O CTSNG projects are facing, this paper studies a coproduction process with LNG and methanol (CTLNG-M). The CTLNG-M process is developed based on a rational distribution study on hydrogen and carbon elements in the processing, which reduce CO2 emission by converting more carbon to chemicals and increase unit product income for a high valued liquefied natural gas (LNG) product. Section 2 gives the description of the new process on what measures have been taken. Section 3 gives the detailed modeling and simulation with respect to key parameters of added units in the CTLNG-M process. In Section 4, a discussion about the carbon utilization efficiency, energy efficiency, energy consumption, and economic performance of the CTLNG-M process is given.

#### **2. LNG and Methanol Coproduction Process**

The syngas from a Lurgi gasifier contains 12% to 18% methane [21]. Because of this high composition of methane, Lurgi gasification technology is usually used in CTSNG projects [22,23]. However, from another point of view, LNG products can also be obtained by separating methane from the syngas through an added cryogenic separation technology. LNG is a relatively high value-added product form of coal-based natural gas, whose price can reflect the supply and demand mechanism. In Figure 3, it can be seen that LNG prices have shown an upward trend from 3206 CNY/tons to 5373 CNY/tons in heating period time (January in each year) since 2017, with an average price of 3122 CNY/tons and a highest price up to 5613 CNY/tons. In addition, LNG can be transported and stored in a more flexible way [24].

Thus, use of a separation unit to remove methane from syngas out of the Lurgi gasifier is taken into consideration, as the remaining syngas can be used for methanol synthesis. In this paper, two units are added, the cryogenic separation unit is placed before a methanol synthesis unit. Thereby, the content of effective gas in methanol synthesis reaction is increased, and the production efficiency also improved.

**Figure 3.** China's liquefied natural gas (LNG) market price recording (January 2017–January 2019).

This paper proposes a coproduction process for matching the product distribution of Lurgi gasification technology. The CTLNG-M process is highlighted through a schematic diagram as shown in Figure 4. In this process, the syngas from the Lurgi gasifier is separated to get LNG product by cryogenic separation unit. The remaining syngas has the H:C ratio close to 2.4. An additional carbon source is needed to decrease that ratio near to 2.1 before methanol synthesis. In that case, the additional carbon resource can be provided by CO2 extracted from the acid gas removal unit. Thereby, the carbon emission of the system is also reduced.

**Figure 4.** Flowchart of coal to liquefied natural gas (LNG) and methanol (CTLNG-M) coproduction process.

There are different gases that are present as impurities in crude syngas. Amongst them, sulfides can cause deactivation of methanol synthesis catalysts, while carbon dioxide can reduce the conversion of methanol synthesis. These impurities, as well as tar, phenol, and ammonia, can be removed by acid gas removal unit [10]. After that and heat recovery, the molecular sieve process is used to further reduce the content of carbon dioxide and methanol to less than 1 ppm and then meet feed requirements of the cryogenic separation [25].

The new process employs nitrogen cycle refrigeration technology to separate the methane [26]. Nitrogen provides most of the cooling capacity through the adiabatic expansion cycle in the turbine expander. A double column cryogenic distillation process is used for separation of syngas and LNG [27], as shown in Figure 5. In the double column process, the washing column and the CH4-CO distillation column are packed columns. The top outlet of the washing tower is syngas with methane content less than 1%. The cold energy is recovered through a heat exchanging system before sending to methanol synthesis and be used for exchange heat from input gas.

**Figure 5.** Process flow diagram of CTLNG-M.

The stream extracted from the bottom of the washing tower mainly consists of methane and carbon monoxide. It is sent to the CH4-CO distillation column for methane separation.

In the CH4-CO distillation column, the condensed liquid stream at the top of the column is partially used as the reflux, and the other part enters the washing column as recycling stream at the top of the column. The main component of the non-condensable gas at the top of this column is CO with the concentration of 70%. In the bottom of tower, a part of LNG returns to the circulation inside the tower for improving product quality with higher CH4 content. The other part is cooled to −163 ◦C through the heat exchanger.

The syngas from the Lurgi gasifier reaches the standard for methanol synthesis through the use of Rectisol and the cryogenic separation unit. After compression, components in syngas react to the product methanol with copper-based catalyst. The main equations are shown as below.

$$\text{CO} + 3\text{H}\_2 = \text{CH}\_3\text{OH} + \text{H}\_2, \,\Delta\text{H} = -90.64 \text{ kJ/mol.} \tag{2}$$

$$\text{CO}\_2 + 3\text{H}\_2 = \text{CH}\_3\text{OH} + \text{H}\_2\text{O}, \Delta\text{H} = -49.47 \text{ kJ/mol.} \tag{3}$$

#### **3. Modeling and Simulation**

As has been mentioned above, there are four main units involved in the CTLNG-M process. Namely coal gasification, acid gas removal, cryogenic separation, and methanol synthesis unit. The detailed simulation of coal gasification and the acid gas removal unit can be found in our group's previous work [28–30]. Consequently, this paper gives modelling and simulation results for two added units, as the cryogenic separation and methanol synthesis unit.

The coal quality parameters are shown in Table 2.

**Table 2.** Proximate and ultimate analyses of Shenmu coal. HHV: high heating value.


#### *3.1. Cryogenic Separation Unit*

The role of the cryogenic separation unit used in the new process is to obtain purified LNG products [31]. As illustrated in Figure 6, the clean syngas is firstly introduced to an absorption tower for H2O and methane removal. After that, the resulted syngas (S1) is cooled into liquid phase. Most of the remaining carbon monoxide component in S2 is sent to the methanol synthesis unit. Stream S4 from the bottom of T2-W mainly consists of CO and CH4, which needs further distillation in T3-D [32]. Stream S5 is then separated in to two parts, one is recycled to the top of T2-W, and another for methanol synthesis. Stream S8 obtained from the bottom of the T3-D is the LNG product.

**Figure 6.** Flow diagram of the cryogenic separation process.

In this paper, the SRK method is used to predict the physical properties of streams. T1-A is modeled by the SEP module, T2-W and T3-F by the RadFrac module and T4-F by the Flash module. The pressure of the top of T2-W is 36 bar, and the theoretical number of tower plates is set to 30, which is twice the minimum theoretical number calculated by Aspen's succinct calculation. The feeding position is stage 14. In the T3-D, its pressure is 34 bar, and drop pressure is 7 bar. In order to maintain the quality of the LNG product, the CH4/(CH4 + C2) ratio of steam gas in the bottom is less than 97.5% (GB/T 19204-2003), which is controlled by adjusting the reflux ratio.

After compression, the nitrogen is liquefied. Through the throttle valve, the high-pressure liquid nitrogen is expanded to a low-pressure state. In this process, the gas absorbs heat from the environment. Therefore, this expansion process provides cooling capacity for the T2-W and the T3-D. The cryogenic system uses two-stage nitrogen circulation expansion refrigeration, which is shown in Figure 7 [33].

**Figure 7.** Flow diagram of the nitrogen circulation refrigeration process.

Simulation results are given in Table 3. Steam S4 consists of 38.8% CH4 and 55.1% CO. Stream S4 is liquefied as the LNG product and outputted from the bottom of T3-D. In the LNG product, the purity of CH4 can approach up to 97.6%, with the impurity content of CO and C2 less than 1%. The separation is modeled and simulated compared with the data in the reference [34], with the error less than 5%.


**Table 3.** Simulation results of the cryogenic separation unit.

#### *3.2. Methanol Synthesis Unit*

The flow chart of the methanol synthesis is shown in Figure 8. Methanol synthesis gas is, at first, mixed with the stream of CO2 from the acid gas removal unit to adjust the H:C ratio. This is because the syngas is rich in hydrogen and additional carbon source for methanol synthesis is needed until that H:C ratio decreased to around 2.1 [35,36].

**Figure 8.** Flow diagram of the methanol synthesis unit.

According to the practical process in Datang Keqi project, the Lurgi low-pressure methanol synthesis method is employed. The methanol syngas enters the reactor with the recycling gas. After heat recovery of the outlet stream, the gas–liquid steam is separated by a separator to recycle the unreacted syngas. The liquid part is the input stream of the methanol rectification unit. A small part of unreacted gas is discharged as purge gas.

In this unit, the PR property method was selected in Aspen for modelling the methanol synthesis unit. The purified syngas is pressurized to 8.2 MPa by the compressor modeled by the COMP module [37]. The synthesis reactor was modeled by a plug flow model of RPlug. Table 4 shows the simulation result of the process.


**Table 4.** Simulation results of the methanol synthesis unit.

The methanol synthesis gas with a H:C ratio of 2.1 is compressed to 82 bar and sent to the reactor [28]. In this syngas, the CO2 volume fraction is less than 3%, which is in line with industrial practice. During the process, 86% unreacted gas is recycled back to the reactor to improve the overall conversion. As a result, the CTLNG-M process produces a total of 6573 kmol/h methanol product.

The key parameters in the CTLNG-M process are listed in Table 5. The carbon element conversion ratio in the coal gasification unit is 99.9%. In the acid gas removal unit, the H2S removal ratio is 99.5%, and the pressure is under 50 bar. In the cryogenic separation unit, the CH4 mole fraction is 97.5%, the recycle ratio and purge ratio in the methanol synthesis unit are set as 0.86 and 0.05 separately.



#### *3.3. Parameters Analysis*

There are two important operation parameters in the new process. One is the temperature and pressure at the bottom of the gas–liquid separation tower in the cryogenic separation unit, and another one is the unreacted gas recycle ratio in the methanol synthesis unit. These two parameters are highly corelated with the composition of syngas and methanol production. These two parameters have significant impact on the composition of syngas and methanol production, which are analyzed in detail as follows.

Figure 9a,b shows the effect of operating temperature and pressure of T4-S on the composition of syngas out from the cryogenic separation unit. For analysis, we fix the recycling ratio to 0.86.

It can be seen from Figure 9a that when the operating pressure of the tower is 31.1 bar, the H:C ratio of syngas decreases from 2.25 to 2.1, and the content of inert gas increases from 0.68% to 0.88% as the temperature rises from −175 to −155 ◦C. When the tower operating pressure is 21.1 bar, as shown in Figure 9b, the H:C ratio is from 2.25 to 1.95, and the content of inert gas increased from 0.68% to 0.95%. The trend of syngas composition changing is similar under different operating pressures.

**Figure 9.** (**a**,**b**) Effect of temperature and pressure on the H:C ratio and inert gas concentration.

However, the H:C ratio of syngas for methanol synthesis should not exceed 2.1 [38]. When the H:C ratio decreases to 2.1 and the pressure is 31.1 bar, the increase of inert gas content is less than that at 21.1 bar. Considering that accumulation of inert gas in syngas will reduce methanol production [17], in this paper the operating pressure is set to 31.1 bar, and the operating temperature is set to −158 ◦C.

We then study the effect of the unreacted gas recycling ratio to the conversion (a) and the compression duty for gas recycling (b), as shown in Figure 10a,b.

The productivity of methanol shows an upward tendency with the recycling ratio increasing, which indicates that a higher carbon utilization efficiency can be achieved by adjusting that ratio. The results confirm a good match with the previous study of Man et al. (2016) [39]. However, with more gas recycled, the units will consume more compression duty, as shown in Figure 10b.

Figure 10a shows that, when the recycle ratio increases from 0.50 to 0.86, the methanol production increases slowly, and when the cycle ratio is more than 0.86, the methanol production increases rapidly and the energy consumption increases rapidly. In order to balance the capacity and energy consumption of the system, in this paper, 0.86 is chosen as the recycle ratio of unreacted methanol gas in this unit.

**Figure 10.** (**a**,**b**) Effect of recycling ratio on methanol productivity and compressor duty.

All results at key points of the CTLNG-M process are shown in Table 6. After dehydration and cooling, the crude syngas of 42,338 kmol/h is sent to the acid gas removal unit. In this unit, the CO2 is removed and the content of CO2 is reduced to 20 ppm and that of the H2S less than 1 ppm. The clean syngas flow (SNYGAS-C) is 28,310 kmol/h with the H:C ratio of 2.47. This stream enters the cryogenic separation unit and is cooled to −168 ◦C. CH4 is separated out of this unit as the LNG product of 4790 kmol/h, denoted by stream S8 in this figure. The yearly output of LNG is 642,000 tons. After separation, CH4 content in the syngas is reduced to 0.6% before entering the methanol reactor. This input stream is mixed in this new process with the CO2 stream from the acid gas removal unit. Finally, the syngas has its H:C ratio at 2.1 and the total yearly methanol output of 1.368 million tons.



#### **4. Discussion**

The CTLNG-M process input consists of 4.656 million tons/a raw coal, which remains the same amount as a benchmark CTSNG process, meanwhile, the outputs consist of 1,367,800 tons/a of methanol and 642,000 tons/a of LNG. The benchmark has the same input amount of coal and outputs consisting of 2 billion Nm3 nature gas only. Based on the simulations, we compare these two processes with respect to energy efficiency, carbon element utilization rate, energy consumption, and economic benefit, as given in Table 7. In the following section, we explain the definition of the indexes and analyze the performances of these two processes.


**Table 7.** System performances parameters of CTLNG-M and CTSNG.

#### *4.1. Energy E*ffi*ciency*

According to the first law of thermodynamics, energy efficiency is defined as the ratio of the energy of effective products (*E*0) to the energy of input raw materials (*Ei*), as given by Equation (4). [40].

$$
\eta = \frac{\sum E\_O}{\sum E\_i} \times 100\%,\tag{4}
$$

where *E*<sup>0</sup> is the product energy (MW) of chemical process and *Ei* is the raw material energy (MW) of chemical process. In this paper, the energy of raw materials and products is calculated by the high heating value (HHV). In the CTLNG-M process, the source of input energy includes raw coal, electricity, and steam. Thereby, the outputs of energy are LNG and methanol. Methanol is a widely-used platform chemical product with a high calorific value of 22.7 GJ/ton. LNG is an energy product and mainly used for urban gas or power generation. It has a high calorific value of 54.6 GJ/ton. Methane is also used as fuel, and the high heating value of the gas conforms to the natural gas GB17820-2012 standard which is 31.4 MJ/m3.

According to Equations (3) and (4), the product energy of the CTLNG-M process is 9158 GJ/h, and the energy efficiency is 53.1%. The product energy of the CTSNG process is 7850 GJ/h, and its energy efficiency is 50.4%. It shows that the new CTLNG-M process has a higher efficiency of 3% than that of the conventional CTSNG process.

#### *4.2. Element Utilization Ratio*

In a coal-based chemical process, the C element in coal is transformed into a chemical product. Thus, it is important to analyze the C resource utilization efficiency to represent the resource utilization. The element C converted into methanol in coal is defined as the effective C, and the element C discharged in the form of CO2 or waste residue is defined as the ineffective C. The ratio of the carbon

mole flow in the product to the mole flow in the raw material is defined as carbon efficiency λ, which can be represented by Equation (5) [41].

$$
\lambda = \frac{\sum F\_O}{\sum F\_i} \times 100\% \tag{5}
$$

where *F*<sup>0</sup> is the mole flow of carbon in methanol and LNG products and *Fi* is the mole flow of carbon in coal.

Figure 11 shows the carbon elemental balance in the new process. It shows that the input molar flow in raw material coal of CTLNG-M system is 27,416 kmol/h. The molar flow of carbon in the crude syngas is 27,141 kmol/h after gasification, and gets 11,749 kmol/h carbon elements when washed with methanol at a low temperature. In the cryogenic separation unit, the molar flow of carbon in the LNG product is 4284 kmol/h. Remaining clean syngas is mixed with the pure CO2 from the acid gas removal unit and the molar flow of carbon in the methanol syngas is 6967 kmol/h. In this case, CO2 through the acid gas removal process, is separated into two parts. Partial CO2 is then removed from gas emission, and recycled in the synthesis process to convert the final methanol product. The methanol syngas re-enters the methanol synthesis unit which contains 6573 kmol/h carbon in the methanol products.

**Figure 11.** Carbon element flow in the process of GTLNG-M.

According to Equation (5), the carbon element efficiency of the new process is 39.6%. This ratio is 5.2% higher than that of the CTSNG process. This is mainly because CO2 emission has been partially converted into product. In the conventional process, all syngas has to be converted to only synthesized natural gas (SNG), which requires the H:C ratio of 3.1 using the element balance equation. This is higher than the ratio in the syngas output from the Lurgi gasification as 2.7. It is necessary for the CTSNG process to use the water–gas shift unit to increase the ratio to 3.1 for methanation reaction. In this course, CO2 emission is increased. However, in the new coproduction process, methanol is present as a suitable product from chemical synthesis through which product methane is separated and cryogenically cooled to directly produce the LNG product. The remaining syngas is only used for methanol synthesis, which requires a lower H:C ratio of 2.1. In this case, the syngas has excessive hydrogen. We then introduce CO2 into the syngas to adjust the ratio. In this study, 209 kmol/h CO2 is converted to methanol.

#### *4.3. Energy Consumption Analysis*

As has been stated in the above discussion, the CTLNG-M process has a higher energy and carbon utilization ratio than CTSNG process. Moreover, considering the new process is under a coproduction design with an added cryogenic separation unit, which is specially needed at low temperature environment and; therefore, consumes more electricity, quantitative analysis for energy use is a necessity.

The energy consumption is defined as utilities consisting of steam cost and electricity cost. According to our calculation, the steam cost in CTLNG-M is 1.34 million tons/a, and the electricity consumption is 110 MW, while the same in the CTSNG process are, respectively, 1.86 million tons/a and

77 MW. For a more convenient comparison, both steam cost and electricity consumption are converted to the same units as MJ/a, as shown in Figure 12.

**Figure 12.** Energy consumption of the CTLNG-M and CTSNG processes.

In Figure 12, the CTLNG-M process consumes 4.3 <sup>×</sup> 109 MJ of steam and 9.5 <sup>×</sup> <sup>10</sup><sup>9</sup> MJ of electricity for a year, and the CTSNG process consumes 5.9 <sup>×</sup> 10<sup>9</sup> MJ of steam and 6.7 <sup>×</sup> 10<sup>9</sup> MJ of electricity. It shows that the coproduction system has a lower steam cost of about 1.6 <sup>×</sup> 109 MJ for per year, which is mainly because of a flexible way to integrate heat exchange when there is not only one route for product processing [42]. However, more electricity is consumed in the new process. It is because the nitrogen circulation refrigeration process needs more power assistance, as modelling data indicates. Since there is no power that can be generated within the system, it takes more capital investment, which needs to be further analyzed.

To summarize, the total energy consumption in general increased by 8.7%. The coproduction process has an advantage on utility usage due to integration of a heat exchanger and flexible distribution flow between different product processing. However, in the specific case of CTLNG-M, a higher electricity consumption is due to compression work in the added cryogenic separation unit. In total, the increased electricity consumption cannot be outweighed by the decrease in the steam cost, and the energy demand gap is 1.2 <sup>×</sup> 10<sup>9</sup> MJ/a, which indicates more investment on various costs in the new coproduction process and a further economic analysis is needed for profitability measurement.

#### *4.4. Economic Analysis*

#### 4.4.1. Total Capital Investment

The total capital investment (TCI) for a given construction project mainly includes fixed capital investment and variable cost. The investment for manufacturing and plant facilities are defined as the fixed capital investment, while those for the plant operation are the working capital [43]. The equipment investment of the system can be calculated by Equations (6) and (7) based on the benchmark investment of the major equipment listed. The total investment can then be derived from the scale factor ( (See Tables A1 and A2 in Appendix.).

$$\mathrm{EI}^{\mathrm{r}} = \sum\_{\mathbf{j}} \boldsymbol{\Theta} \cdot \mathrm{EI}\_{\mathbf{j}}^{\mathrm{r}} \cdot (\frac{\mathbf{S}\_{\mathbf{j}}}{\mathbf{S}\_{\mathbf{j}}^{\mathrm{r}}})^{\mathrm{sf}},\tag{6}$$

$$\text{TCI} = \text{EI} \{ 1 + \sum\_{i=1}^{n} RF\_i \} \tag{7}$$

where EI is the equipment investment (CNY), θ is the localization factor, EI<sup>r</sup> <sup>j</sup> is the benchmark equipment investment of the j unit, Sj is the scale of the j unit, S <sup>r</sup> <sup>j</sup> is the base scale of the j unit, TCI is the total investment, and RFi is the proportionality factor of the investment composition *i*.

As shown in Figure 13, based on the same input of raw coal, the total investment of the CTSNG project is 16.62 billion CNY, and the CTLNG-M project is 13.66 billion CNY, which is 17.8% lower than CTSNG. This is because the new process eliminates the water–gas shift unit compared to the single-production coal gasification process, so that the carbon emission is less and the amount of gas processed by the acid gas removal unit is decreasing compared to the original CTSNG process. The corresponding investment is also reduced. At the same time, the CTLNG-M process uses a nitrogen expansion refrigeration process, with mature technology and low investment. Therefore, the cryogenic unit equipment and related investment are relatively low, and the total amount of process investment is correspondingly reduced, which is more suitable for CTSNG projects and has economic advantages.

**Figure 13.** Capital investment of CTLNG-M and CTSNG.

#### 4.4.2. Internal Rate of Return

Internal rate of return (IRR) is another important index for evaluation of economic performance, which takes into account the net present value and the service life of processing route into account [44]. A dynamic evaluation method is taken in this paper, the calculation is as follows.

$$\text{NPV} = \sum\_{t=0}^{\text{m}} \frac{(CI - CO)\_t}{(1 + i)^t} \,\text{},\tag{8}$$

$$\text{NPV} = \sum\_{t=0}^{m} \frac{(CI - CO)\_t}{\left(1 + IRR\right)^t} = 0,\tag{9}$$

where *CI* is the net cash inflow in the *t* year, *CO* is the net cash outflow; *m* is the project's life time; *i* is the benchmark rate of return. NPV stands for net present value (NPV), which refers to the net cash flow generated annually by a technical solution throughout its life cycle. The net cash flow generated each year is converted to the present value at the base time by a specified base discount rate *i*0.

Inner rate of return (IRR) can usually be calculated by interpolation method. It represents the discount rate when the cumulative present value of the net cash flow of the project is equal to zero in the whole calculation period. IRR is a dynamic index to evaluate the economic feasibility of new projects. It is usually compared with the base rate of return to determine whether the new chemical process is feasible. In this paper, *i* is set to 12% [45]. If the IRR is larger than the base rate of return *i*, the process is economically feasible and achieves the lowest level of return on investment. In addition, with the increase of internal rate of return, the obtained benefit of the process will also increase.

The IRRs of the CTLNG-M and CTSNG processes are compared in Figure 14, which are higher than the industrial criterion of 12%. Specifically, the IRR of the CTSNG process is 13%, which is slightly higher than 12%, and accords with the current status of the CTSNG project. However, the IRR of the CTLNG-M process is 19%, which increased by 6%, so this process has higher profit.

**Figure 14.** Inner rate of return (IRR) of the CTLNG-M and CTSNG processes.

#### **5. Conclusions**

This paper proposes a system of coproduction for LNG and methanol. The aim was to find improvements to the low-earning CTSNG process using the same raw material but producing a low-margin, single SNG product. In the new coproduction process, there are two innovative aspects. On the one hand, the syngas is firstly separated to the LNG product and the lean-methane syngas is then used for methanol synthesis. To realize this improvement, a cryogenic separation unit is added. Besides, the syngas with little CH4 has a higher hydrogen component than that for methanol synthesis. Thereby, CO2 is used to supply an additional carbon element to the methanol synthesis. On the other hand, the methanation unit is removed, while the process still outputs a product of the high-valued form of methane as the LNG. In the case study, we modeled and simulated the key units of the CTLNG-M process with 642,000 tons/a LNG and 1.368 million tons/a methanol product, compared to the CTSNG process with the same coal processing coal capacity and 2 billion NM3/a SNG. In element efficiency analysis, the carbon efficiency of the new process increases from 34.7% to 39.6%, with corresponding decrease of carbon emission by 130,000 tons per year. Because of the additional energy consumption for gas compression, the energy efficiency of the new process is at the same level with the CTSNG process.

In economic analysis, the IRR of the CTSNG process is 13% while the IRR of the CTLNG-M process is 19%. The new process brings much higher economic benefits. This is because the new process produces a higher valued product and saves the carbon resource during methanol synthesis. Moreover, the new process has 17.8% reduction of investment compared to the CTSNG process. Thus, this is a promising solution for coal chemical processes based on Lugri gasification technology, with more economic benefit and less investment.

**Author Contributions:** J.G. and S.Y. contributed to the conception of the study. J.G. and S.Y. contributed significantly to analysis and manuscript preparation; J.G. performed the data analyses and wrote the manuscript; A.K. helped perform the analysis with constructive discussions.

**Funding:** We would like to express our appreciation to the National Natural Science Foundation of China (funder ID: http://dx.doi.org/10.13039/501100001809, Grant No. 21736004).

**Conflicts of Interest:** The authors declare no conflicts of interest.

#### **Appendix A**


**Table A1.** Summary of investment data for main equipment components.



#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Statistical Process Monitoring of the Tennessee Eastman Process Using Parallel Autoassociative Neural Networks and a Large Dataset**

#### **Seongmin Heo and Jay H. Lee \***

Department of Chemical and Biomolecular Engineering, Korea Advanced Institute of Science and Technology (KAIST), 291 Daehak-ro, Yuseong-gu, Daejeon 34141, Korea

**\*** Correspondence: jayhlee@kaist.ac.kr

Received: 3 May 2019; Accepted: 21 June 2019; Published: 1 July 2019

**Abstract:** In this article, the statistical process monitoring problem of the Tennessee Eastman process is considered using deep learning techniques. This work is motivated by three limitations of the existing works for such problem. First, although deep learning has been used for process monitoring extensively, in the majority of the existing works, the neural networks were trained in a supervised manner assuming that the normal/fault labels were available. However, this is not always the case in real applications. Thus, in this work, autoassociative neural networks are used, which are trained in an unsupervised fashion. Another limitation is that the typical dataset used for the monitoring of the Tennessee Eastman process is comprised of just a small number of data samples, which can be highly limiting for deep learning. The dataset used in this work is 500-times larger than the typically-used dataset and is large enough for deep learning. Lastly, an alternative neural network architecture, which is called parallel autoassociative neural networks, is proposed to decouple the training of different principal components. The proposed architecture is expected to address the co-adaptation issue of the fully-connected autoassociative neural networks. An extensive case study is designed and performed to evaluate the effects of the following neural network settings: neural network size, type of regularization, training objective function, and training epoch. The results are compared with those obtained using linear principal component analysis, and the advantages and limitations of the parallel autoassociative neural networks are illustrated.

**Keywords:** process monitoring; nonlinear principal component analysis; parallel neural networks; autoassociative neural network; big data

#### **1. Introduction**

Statistical process monitoring is one of the most intensely-studied problems for the modern process industry. With the rising need for sustainable operation, it has been attracting extensive research effort in the last few decades [1,2]. The key step of statistical process monitoring is to define normal operating regions by applying statistical techniques to data samples obtained from the process system. Typical examples of such techniques include principal component analysis (PCA) [3–6], partial least squares [7–9], independent component analysis [10,11], and support vector machine [12,13]. Any data sample that does not lie in the normal operating region is then classified as a fault, and its root cause needs to be identified through fault diagnosis.

Recently, deep learning and neural networks have been widely used for the purpose of statistical process monitoring, where both supervised and unsupervised learning algorithms have been implemented. In designing process monitoring systems in a supervised manner, various types of neural networks have been used, such as feedforward neural networks [14], deep belief networks [15], convolutional neural networks [16], and recurrent neural networks [17]. In the case of unsupervised

learning, the autoassociative neural network (also known as an autoencoder), which has been proposed as a nonlinear generalization of PCA [18], is typically used [19–21]. While the traditional statistical approaches typically rely only on the normal operating data to develop the process monitoring systems, most of the deep learning-based process monitoring studies have adopted supervised learning approaches. However, in the real industrial processes, it is difficult to obtain a large number of data samples for different fault types, which can be used for the training of deep neural networks. Thus, it is important to examine rigorously the potential of autoassociative neural networks as a basis for the design of process monitoring systems.

In the process systems area, the Tennessee Eastman (TE) process, a benchmark chemical process introduced by Downs and Vogel [22], has been a popular test bed for process monitoring techniques. There already exist a few studies where the process monitoring systems for this process are designed on the basis of autoassociative neural networks [23–25]. However, considering the complexity of the neural network training, these studies have two limitations. First, a rigorous case study has not been performed to evaluate the effects of different neural network settings, such as neural network hyperparameters and training objective functions, which can have great impact on the performance of the process monitoring systems. Furthermore, a few thousand normal training samples were used to train neural networks with much more parameters, ranging from a hundred thousand to a million parameters. A larger dataset is required to investigate the effectiveness of unsupervised deep learning for the statistical process monitoring. In addition to the above limitations, there is another issue that is directly related to the structure of autoassociative neural networks. It has been reported that there is a high chance for the principal components, which are extracted using autoassociative neural networks, to be redundant due to the co-adaptation in the early phase of neural network training [18]. The objective of this work is to address these limitations.

The rest of the article is organized as follows. First, the concept of linear PCA is briefly explained, and how it can be used for the statistical process monitoring is discussed. Then, the information on autoassociative neural networks is provided, and the parallel autoassociative neural network architecture, which was proposed in our previous work [26] to alleviate the co-adaptation issue mentioned above, is described. This is followed by the description of the statistical process monitoring procedure using autoassociative neural networks. Finally, a comprehensive case study is designed and performed to evaluate the effects of different neural network settings on the process monitoring performance. The dataset used in this study has 250,000 normal training samples, which is much larger than the ones considered in the previous studies.

#### **2. Principal Component Analysis and Statistical Process Monitoring**

#### *2.1. Linear Principal Component Analysis*

Let us first briefly review the concept of linear PCA. PCA is a statistical technique that decorrelates the original variables, resulting in a set of uncorrelated variables called principal components. Let *x* be a sample vector of *m* variables and *X* be a data matrix whose rows represent *n* sample vectors. Assuming that each column of *X* has zero mean and unit variance, the singular value decomposition can be applied to the sample covariance matrix:

$$\frac{1}{m-1}X^TX = P\Lambda P^T\tag{1}$$

where Λ is a diagonal matrix that contains the eigenvalues of the sample covariance matrix on its main diagonal and *P* denotes an orthogonal matrix whose columns are the eigenvectors of the sample covariance matrix. In the context of PCA, *P* is called the loading matrix since its column vectors can be used to extract principal components from the original data as follows:

$$T = XP\tag{2}$$

where *T* represents the score matrix whose elements are the principal component values.

Let λ*<sup>i</sup>* be the diagonal element in the *i* th row of Λ, and let us assume that Λ is arranged in a descending order (i.e., λ<sup>1</sup> ≥ λ<sup>2</sup> ≥···≥ λ*m*) so that the *j* th column of *P* corresponds to the direction with the *j* th largest variance in the principal component space. Then, we can partition the loading matrix into two blocks as below:

$$P = \begin{bmatrix} P\_{PC} & P\_R \end{bmatrix} \tag{3}$$

where *PPC* and *PR* contain the first *l* columns and the remaining columns of *P*, respectively. These two submatrices can be respectively used to map the original data onto the lower dimensional principal component space and residual space (of dimension *l* and dimension *n-l*, respectively):

$$\begin{aligned} T\_{\text{PC}} &= X P\_{\text{PC}} \\ T\_R &= X P\_R \end{aligned} \tag{4}$$

where *TPC* and *TR* are the first *l* columns and the remaining columns of *T*, respectively. In what follows, we explain how linear PCA can be used for statistical process monitoring.

#### *2.2. Statistical Process Monitoring Using Linear PCA*

In the PCA-based process monitoring, a new data sample is first projected onto the lower dimensional principal component space and the residual space [27]. Then, it is evaluated whether the new sample lies in the normal operating range in both spaces. Hotelling's *T2* and *Q* (or squared prediction error) statistics are typically used to define the normal operating range in the principal component space and the residual space, respectively, for such evaluation. These statistics can be computed by the following equations:

$$\begin{array}{c} T^2 = \mathbf{x} P\_{\mathbf{PC}} \boldsymbol{\Lambda}\_{\mathbf{PC}}^{-1} P\_{\mathbf{PC}}^T \mathbf{x}^T \\ Q = \mathbf{x} P\_{\mathbf{R}} P\_{\mathbf{R}}^T \mathbf{x}^T \end{array} \tag{5}$$

where Λ*PC* represents a diagonal matrix formed by the first *l* rows and columns of Λ.

The upper control limit for the *T2* statistic is given as [28]:

$$T\_{\alpha}^{2} = \frac{l(n^{2} - l)}{n(n - l)} F\_{\alpha}(l, n - l) \tag{6}$$

where *F*α(*l*, *n* − *l*) represents the α percentile of the *F*-distribution with *l* and *n-l* degrees of freedom. The *Q* statistic has the upper limit of the following form [29]:

$$Q\_a = \lg \chi\_a^2(h) \tag{7}$$

where:

$$\begin{array}{lcl} \mathbf{g} = \theta\_2 / \theta\_1 \\ \mathbf{h} = \theta\_1^2 / \theta\_2 \\ \theta\_i = \sum\_{k=l+1}^{m} \lambda\_{k'}^i & i = 1, 2 \end{array} \tag{8}$$

and χ<sup>2</sup> <sup>α</sup>(*h*) is the α percentile of the χ*2*-distribution with *h* degrees of freedom.

If a very large number of data samples is available, the mean and covariance estimated from the data will be very close to the true values of the underlying probability distribution. In this case, the upper control limit for *T2* statistic takes the following form [30]:

$$T\_a^2 = \chi\_a^2(l) \tag{9}$$

while the upper control limit for the *Q* statistic can be approximated by the following equation [29]:

$$\begin{array}{l} Q\_{\alpha} = \text{g}\chi^{2}\_{\alpha}(h) \\ g = \frac{\sigma^{2}\_{R}}{2\mu\_{R}} \\ h = \frac{2\mu\_{R}^{2}}{\sigma\_{R}^{2}} \end{array} \tag{10}$$

where μ*<sup>R</sup>* and σ*<sup>R</sup>* are the mean and standard deviation of the squared prediction errors (i.e., *Q* values) obtained from the training dataset.

If any of the above control limits is violated, the new sample is classified as a fault. Once a fault is detected, its root cause needs to be identified. The contribution plot is typically used to this end, where the contribution of each variable to the *T2* or *Q* statistic is calculated [31–33].

It is also important to select a proper number of principal components to be retained in the PCA model. For such selection, various criteria are available in the literature including cumulative percent variance [34], residual percent variance [35], parallel analysis [36], and cross-validation [37]. For a detailed discussion on this subject, the readers are referred to the work by Valle et al. [38].

#### **3. Statistical Process Monitoring Using Autoassociative Neural Networks**

#### *3.1. Nonlinear Principal Component Analysis Using Autoassociative Neural Networks*

We now describe neural network-based nonlinear principal component analysis (NLPCA). Kramer [18] proposed to use a special type of neural network, called the autoassociative neural network, for NLPCA. As shown in Figure 1, an autoassociative neural network consists of five layers: input, mapping, bottleneck, demapping, and output layers. Its goal is to learn the identity mapping function to reconstruct its input data at the output layer. The problem of learning identity mapping becomes non-trivial if the dimension of the bottleneck layer, *f*, is smaller than that of the original data, *m*.

**Figure 1.** Network architecture of the autoassociative neural network.

The first three layers of autoassociative neural network approximate the mapping functions, which project the original data onto the lower dimensional principal component space, while the last

two layers approximate the demapping functions, which bring back the projected data to the original data space. The mathematical model of the autoassociative neural network has the following form:

$$\begin{aligned} y\_m &= a(x\mathcal{W}\_1 + b\_1) \\ t &= y\_m \mathcal{W}\_2 + b\_2 \\ y\_d &= a(t\mathcal{W}\_3 + b\_3) \\ \mathfrak{k} &= y\_d \mathcal{W}\_4 + b\_4 \end{aligned} \tag{11}$$

where *x*, *ym*, *t*, *yd*, and *x*ˆ represent the vectors of input, mapping, bottleneck, demapping, and output layers, respectively. *W* and *b* are weight matrices and bias vectors, respectively. The dimensions of all the matrices and vectors are summarized in Figure 1. *a* denotes the nonlinear activation function. The objective of autoassociative neural network training is to find optimal parameter values (i.e., optimal values of *W* and *b*) that minimize the difference between the input and the output, i.e.:

$$E = \frac{\sum\_{i=1}^{n} \sum\_{j=1}^{m} \left(\mathbf{x}\_{ij} - \mathbf{x}\_{ij}\right)^2}{nm} \tag{12}$$

which is also called the reconstruction error.

#### *3.2. Alternative Neural Network Architecture: Parallel Autoassociative Neural Networks*

It was pointed out by Kramer [18] that principal components extracted from an autoassociative neural network can be redundant, as multiple principal components are aligned together in the early stage of network training. To this end, in our previous work [26], we proposed an alternative neural network architecture to address this limitation, which decouples the training of different principal components. Such decoupling can be achieved by alternating the network architecture of the autoassociative neural network as shown in Figure 2.

**Figure 2.** Alternative neural network architecture for parallel extraction of principal components.

In this network architecture, all the hidden layers are decomposed into *f* sub-layers to form *f* decoupled parallel subnetworks. Each subnetwork extracts one principal component directly from the input and reconstructs the pattern that it captures. Then, the outputs from all the subnetworks are added up to reconstruct the input. The mathematical model of this network architecture has the same form as the one in Equation (11) with the structural changes to the weight matrices and bias vectors as below:

$$\begin{aligned} W\_1 &= \begin{bmatrix} W\_{11} & W\_{12} & \cdots & W\_{1f} \\ W\_{21} & 0 & \cdots & 0 \\ 0 & W\_{22} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & W\_{2f} \\ 0 & 0 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & W\_{3f} \end{bmatrix} \\ W\_4 &= \begin{bmatrix} W\_{41} \\ W\_{42} \\ \vdots \\ W\_{4f} \end{bmatrix} \\ b\_1 &= \begin{bmatrix} b\_{11} & b\_{12} & \cdots & b\_{1f} \\ & \vdots \\ b\_{21} & b\_{22} & \cdots & b\_{2f} \\ \vdots \\ b\_{31} & b\_{32} & \cdots & b\_{3f} \end{bmatrix} \\ b\_4 &= b\_{41} + b\_{42} + \cdots + b\_{4f} \end{aligned} \tag{13}$$

The readers are referred to [26] for more detailed information on the parallel autoassociative neural networks (e.g., the systematic approach for the network decoupling and the advantages of the decoupled parallel neural networks). The proposed network architecture has two potential advantages over the existing one, which are relevant to the statistical process monitoring. First, due to the decoupling, the proposed network architecture is expected to extract more independent (i.e., less correlated) principal components and result in smaller reconstruction errors compared to the existing architecture. If we can achieve smaller reconstruction error using the same number of principal components, it can imply that the more essential information of the original data is captured. Thus, there is a potential for small reconstruction error to translate into high process monitoring performance. The other advantage is that the proposed network architecture requires much fewer parameters compared to the existing architecture, given that the networks have the same size (i.e., the same number of hidden layers and nodes). As a result, the proposed architecture is expected to be more robust to network overfitting, which can lead to more consistent process monitoring performance. Furthermore, the proposed architecture is more suitable for online implementation since it can compute the values of the *T<sup>2</sup>* and *Q* statistics more quickly than the existing architecture.

#### *3.3. Objective Functions for Autoassociative Neural Network Training*

Besides the reconstruction error in Equation (12), there exist several objective functions available for the autoassociative neural network training. Here, we provide a brief description of two alternative objective functions: hierarchical error and denoising criterion. Hierarchical error was proposed by Scholz and Vigário [39] to develop a hierarchy (i.e., relative importance) among the nonlinear principal components, which does not generally exist for the principal components obtained by using the reconstruction error as the objective function. In linear PCA, it can be shown that the maximization of

the principal component variance is equivalent to the minimization of the residual variance. Motivated by this, the following hierarchical reconstruction error can be defined:

$$E\_H = \sum\_{k=1}^f \alpha\_k E\_k \tag{14}$$

where *Ek* represents the reconstruction error calculated by using only the first *k* nodes in the bottleneck layer and α*<sup>k</sup>* is a hyperparameter that balances the trade-off among the different error terms. The problem of selecting the optimal values of α*<sup>k</sup>* can be computationally expensive, especially in the case of a large bottleneck layer (i.e., large number of principal components). It was illustrated that setting the values of α*<sup>k</sup>* to one can robustly balance the trade-off among the different error terms [39].

The denoising criterion was proposed by Vincent et al. [40] to extract more robust principal components. To apply the denoising criterion, the corrupted input '*x* is generated by adding a noise, such as Gaussian noise and masking noise, to the original input *x*. Then, the autoassociative neural network is trained such that it can recover the original input from the corrupted input. It was shown that, using the denoising criterion, autoassociative neural networks were able to learn a lower dimensional manifold that captures more essential patterns in the original data. The three objective functions are schematically summarized in Figure 3.

**Figure 3.** Schematic representation of different objective functions: (**a**) reconstruction error; (**b**) hierarchical error; (**c**) denoising criterion.

#### *3.4. Statistical Process Monitoring Using NLPCA*

We can design a similar procedure for statistical process monitoring using autoassociative neural networks. First, the original data matrix, which contains only the normal operating data, is partitioned into two disjoint sets, one for the training and the other for the testing of neural networks. Then, an autoassociative neural network is trained using the training dataset. Once the network training is complete, a new data sample is provided to the trained autoassociative neural network to compute principal components and residuals. The *T2* and *Q* statistics can then be calculated as follows:

$$\begin{aligned} T^2 &= \sum\_{k=1}^f \frac{t\_k^2}{\sigma\_k^2} \\ Q &= (\mathbf{x} - \hat{\mathbf{x}})(\mathbf{x} - \hat{\mathbf{x}})^T \end{aligned} \tag{15}$$

where *tk* is the value of the *k*th principal component for the new data sample and σ*<sup>k</sup>* represents the standard deviation of the *k*th principal component calculated from the training dataset.

Note that the upper control limits presented in the previous section are obtained by assuming that the data follow a multivariate normal distribution. In the case of linear PCA, if the original data are normal random vectors, the principal components and residuals also have multivariate normal distributions. Thus, the limits in Equations (6)–(10) can be directly applied to the statistical process monitoring using linear PCA. However, in the case of NLPCA, it is not guaranteed that the principal components follow a multivariate normal distribution since they are obtained by nonlinear transformations. Therefore, in this work, we take an alternative approach, where the upper control limits for two statistics are directly calculated from the data without assuming a specific type of probability distribution, given that a large dataset is available. For example, with 100 normal training data samples, the second largest *T<sup>2</sup>* (or *Q*) value is selected to be the upper control limit to achieve the false alarm rate of 0.01.

#### **4. Process Monitoring of the Tennessee Eastman Process**

Let us now evaluate the performance of the NLPCA-based statistical process monitoring with the Tennessee Eastman (TE) process as an illustrative example. The TE process is a benchmark chemical process [22], which involves five major process units (reactor, condenser, compressor, separator, and stripper) and eight chemical compounds (from A–H). A data sample from this process is a vector of 52 variables, and there are 21 programmed fault types. In this work, Faults 3, 9, and 15 are not considered since they are known to be difficult to detect due to no observable change in the data statistics [41]. The large dataset provided by Rieth et al. [42] is used in this study, and the data structure is summarized in Table 1. Note that Fault 21 is not included in this dataset, and thus not considered in this study. This dataset includes data samples from 500 simulation runs as the training data of each operating mode (normal and 20 fault types) and from another 500 simulation runs as the test data of each operating mode. From each simulation run, 500 data samples were obtained for training, while 960 data samples were recorded for testing. Different types of faults were introduced to the process after Sample Numbers 20 and 160 for the fault training and fault testing, respectively.


**Table 1.** Number of samples in each data subset.

All the neural networks were trained for 1000 training epochs with the learning rate of 0.001. The rectified linear unit (ReLU) was used as the nonlinear activation function, which is defined as max(0,x). The ADAM optimizer [43] and the Xavier initialization [44] were used for the network training. The results reported here are the average values of 10 simulation runs.

In what follows, we first check the validity of the upper control limits estimated using the data only. Then, the performance of the process monitoring using NLPCA is evaluated by analyzing the effects of various neural network settings. Finally, the performance of the NLPCA-based process monitoring is compared with the linear-PCA-based process monitoring.

#### *4.1. Upper Control Limit Estimation*

Let us first compare the upper control limits calculated from the *F*- and χ*2*-distributions and from the data distribution. The objective of this analysis is to show that the dataset used in this study is large enough so that the upper control limits can be well approximated from the data. Another dataset, which contains only 500 normal training samples, is also used for illustration purposes and

was obtained from http://web.mit.edu/braatzgroup/links.html. This dataset and the large dataset from Rieth et al. [42] will be denoted as the S (small) dataset and L (large) dataset, respectively.

Equations (9) and (10) were used to calculate the upper control limits for the L dataset on the basis of the *F*- and χ*2*-distributions, while those for the S dataset were computed using Equations (6)–(8). Linear PCA was used to calculate the upper control limits in the principal component space and the residual space with α = 0.99, and the results are tabulated in Table 2. Note that, for the L dataset, the control limits obtained directly from the data had almost the same values as the ones from the *F*- and χ*2*-distributions, while large deviations were observed for the S dataset. Thus, in the subsequent analyses, the control limits obtained directly from the data will be used for both linear PCA and nonlinear PCA.


**Table 2.** Upper control limits calculated from the probability distributions and the data.

#### *4.2. Neural Network Hyperparameters*

In this work, two neural network architectures were used to evaluate the NLPCA-based process monitoring. The NLPCA methods utilizing the networks shown in Figures 1 and 2 will be referred to as sm-NLPCA (simultaneous NLPCA) and p-NLPCA (parallel NLPCA), respectively. The performance of the process monitoring was evaluated by two indices, fault detection rate (FDR) and false alarm rate (FAR). The effects of the following hyperparameters were analyzed:


We designed four different types of neural networks to evaluate the effects of the hyperparameters listed above. Types 1 and 2 had five layers (three hidden layers), while seven layers (five hidden layers) were used for Types 3 and 4. In Types 1 and 3, the numbers of mapping/demapping nodes were fixed at specific values, while they were proportional to the number of principal components to be extracted for Types 2 and 4. The number of parameters for different network types are summarized in Table 3. The numbers for the network structures represent the number of nodes in each layer starting from the input layer. Note that, for the same network type, p-NLPCA always had fewer parameters compared to sm-NLPCA due to the network decoupling.


**Table 3.** Number of parameters for different neural network types.

Tables 4 and 5 show the process monitoring results for Types 1 and 2 and Types 3 and 4, respectively. Note that, in this analysis, only the average value of FDR (over all the fault types) is reported for brevity.


**Table 4.** Process monitoring results with varying the neural network size (with five layers).




The main trends to note are:


Regarding the last point, in the case of the demapping functions, the input had a lower dimension than the output, which made the problem of approximating demapping functions ill-posed. Thus, it can be speculated that the network overfitting mainly occurred during the reconstruction of the data (i.e., demapping functions were overfitted), leaving the results in the principal component space unaffected by the network overfitting. In addition to this, it was observed that, by including more nodes in the mapping/demapping layers, the average standard deviation of the principal components was increased by a factor of 2~8. This implies that, in Network Types 2 and 4, the normal operating region in the principal component space was more "loosely" defined (i.e., the normal data cluster had a larger volume) compared to Network Types 1 and 3, which can make the problem of approximating the demapping functions more ill-posed. It can also be a reason why the FAR in the principal component space was consistently low regardless of the network size and the degree of network overfitting.

Table 6 shows the FDR values obtained by adjusting the upper control limits such that the FAR became 0.01 for the normal test data. The following can be clearly seen:

• sm-NLPCA performed better than p-NLPCA in the principal component space, while p-NLPCA was better in the residual space.


**Table 6.** Fault detection rates with adjusted upper control limits (FAR = 0.01 for normal test data).

By comparing the results from Network Types 1 and 3, adding additional hidden layers was shown to improve the FDR in the principal component space for sm-NLPCA and the FDR in the residual space for p-NLPCA. However, the effects of such addition cannot be evaluated clearly for Network Types 2 and 4. Thus, in what follows, we apply the neural network regularization techniques to Network Types 2 and 4 and evaluate the effects of such techniques on the performance of NLPCA-based process monitoring.

#### *4.3. Neural Network Regularization*

In this analysis, we consider three different types of neural network regularization: dropout and *L1* and *L2* regularizations. Dropout is a neural network regularization technique, where randomly-selected nodes and their connections are dropped during the network training to prevent co-adaptation [47]. *L1* and *L2* regularizations prevent network overfitting by putting constraints on the *L1* norm and *L2* norm of the weight matrices, respectively.

The process monitoring results for Network Types 2 and 4 are tabulated in Tables 7 and 8, respectively. In the case of Network Type 2, dropout did not address the problem of overfitting for p-NLPCA, and therefore, the results obtained using dropout are not presented here.


**Table 7.** Process monitoring results with neural network regularization (Network Type 2).

**Table 8.** Process monitoring results with neural network regularization (Network Type 4).



• In most cases, neural network regularization degraded the process monitoring performance in the principal component space with p-NLPCA of Network Type 4 regularized by dropout being the only exception.

• On the other hand, it dramatically reduced the FAR values in the residual space for sm-NLPCA, while such reduction was not significant for p-NLPCA (the FAR in the residual space even increased for Network Type 2). This indicates that overfitting was a problem for the residual space detection using sm-NLPCA.

Table 9 shows the process monitoring results obtained by adjusting the upper control limits as mentioned in the previous analysis. Overall, p-NLPCA performed better than sm-NLPCA, while sm-NLPCA was better than p-NLPCA in the principal component space when Network Type 4 was used. By comparing the results provided in Tables 6 and 9, it can be seen that having more nodes in the mapping/demapping layers is only beneficial to the process monitoring in the principal component space. Putting more nodes implies the increased complexity of the functions approximated by neural networks. Thus, it makes the problem of approximating demapping functions more ill-posed and has the potential to be detrimental to the performance of process monitoring in the residual space. Although neural network regularization techniques can reduce the stiffness and complexity of the functions approximated by neural networks [48], they seem to be unable to define better boundaries for one-class classifiers.

**Table 9.** Process monitoring results with network regularization and adjusted upper control limits (FAR = 0.01 for the normal test data).


#### *4.4. Network Training Objective Function*

Let us now evaluate the effects of different objective functions on the process monitoring performance. For illustration purposes, only Network Type 4 was considered in this analysis. All the values of α*<sup>k</sup>* were set to one for the hierarchical error, and the corrupted input was generated by using a Gaussian noise of zero mean and 0.1 standard deviation for the denoising criterion. *L2* regularization and *L1* regularization were used to prevent overfitting for sm-NLPCA and p-NLPCA, respectively. Tables 10 and 11 summarize the process monitoring results obtained by using the autoassociative neural networks trained with different objective functions.

The following are the major trends to note:


**Table 10.** Process monitoring results with different neural network training objective functions.


**Table 11.** Process monitoring results with different neural network training objective functions and adjusted upper control limits (FAR = 0.01 for normal test data).


#### *4.5. Neural Network Training Epochs*

In this analysis, the effects of neural network training epochs are analyzed. To this end, Network Type 4 with 15 principal components was trained, and the neural network parameters were saved at every 10 epochs. Figure 4 shows how different values evolve as the neural networks are trained. For the reference case, the neural networks were trained without any regularization and with the reconstruction error as the objective function. It can be clearly seen that the network overfitting (which is captured by the difference between the solid and dashed black lines) resulted in high FAR values in the residual space (which is captured by the difference between the solid and dashed red lines), while

it did not affect the FAR values in the principal component space (which is captured by the difference between the solid and dashed blue lines).

**Figure 4.** Change in the process monitoring results with respect to the training epoch.

Note the following:


Thus, during the network training, it was required to monitor both reconstruction error and process monitoring performance indices, and early stopping needed to be applied as necessary to ensure high monitoring performances. Nonetheless, from the above observations, it can be concluded that the objective functions available in the literature, which focus on the reconstruction ability of the autoassociative neural networks, may not be most suitable for the design of one-class classifiers. This necessitates the development of alternative training algorithms of autoassociative neural networks to improve the performance of neural-network-based one-class classifiers.

#### *4.6. Comparison with Linear-PCA-Based Process Monitoring*

Let us finally compare the NLPCA-based process monitoring with the linear-PCA-based one. For the linear-PCA-based process monitoring, based on the parallel analysis [35], the number of principal components to be retained was selected as 15. The same number of principal components was used in the NLPCA-based process monitoring. For sm-NLPCA, the following setting was used: Network Type 2, no regularization, reconstruction error as the objective function. For p-NLPCA, the following setting was used: Network Type 3, no regularization, reconstruction error as the objective function.

Table 12 tabulates the process monitoring results obtained using three different PCA methods. Note that the upper control limits were adjusted to have FAR values of 0.01. The main trends observed were:


Let us consider two cases that illustrate the advantages of p-NLPCA over the linear PCA as the basis for the process monitoring system design. Figure 5 shows the *Q* statistic values (black solid line) for one complete simulation run with Fault 5, along with the upper control limit (dashed red line). Although both linear PCA and p-NLPCA detected the fault very quickly (fault introduced after Sample Number 160 and detected at Sample Number 162), in the case of the linear PCA, the *Q* statistic value dropped below the upper control limit after around Sample Number 400. The *Q* statistic value calculated from p-NLPCA did not decrease much, indicating that the fault was not yet removed from the system.


**Table 12.** Detailed process monitoring results obtained using linear PCA and two NLPCA methods (FAR = 0.01 for normal test data).

The contribution plots of *Q* statistic for Fault 1, which involves a step change in the flowrate of the A feed stream, are provided in Figure 6. In the case of the linear PCA, the variables with the highest contribution to the *Q* statistic were Variables 4 and 6, which are the flowrates of Stream 4 (which contains both A and C) and the reactor feed rate, respectively. Note that although these variables were also highly affected by the fault, they were not the root cause of the fault. On the other hand, in the case of p-NLPCA, the variables with the highest contribution to the *Q* statistic were Variables 1 and 44, which both represent the flowrate of the A feed stream, the root cause of the fault. Thus, it can be concluded that the process monitoring using p-NLPCA showed some potential to perform better at identifying the root cause of the fault introduced to the system.

**Figure 5.** Process monitoring results in the residual space for Fault 5: (**a**) linear PCA; (**b**) p-NLPCA.

**Figure 6.** Contribution plots of *Q* statistics for Fault 1.

#### **5. Conclusions**

The statistical process monitoring problem of the Tennessee Eastman process was considered in this study using autoassociative neural networks to define normal operating regions. Using the large dataset allowed us to estimate the upper control limits for the process monitoring directly from the data distribution and to train relatively large neural networks without overfitting. It was shown that the process monitoring performance was very sensitive to the neural network settings such as neural network size and neural network regularization. p-NLPCA was shown to be more effective for the process monitoring than the linear PCA and sm-NLPCA in the residual space, while its performance was worse than the others in the principal component space. p-NLPCA also showed the potential of better fault diagnosis capability than the linear PCA, locating the root cause more correctly for some fault types.

There still exist several issues that need to be addressed to make autoassociative neural networks more attractive as a tool for statistical process monitoring. First, a systematic procedure needs to be developed to provide a guideline for the design of optimal autoassociative neural networks to be used for the statistical process monitoring. Furthermore, a new neural network training algorithm may be required to extract principal components that are more relevant to the process monitoring tasks. Finally, the compatibility of different techniques to define the upper control limits, other than the *T<sup>2</sup>* and *Q* statistics, needs to be extensively evaluated.

**Author Contributions:** Conceptualization, S.H. and J.H.L.; formal analysis, S.H.; funding acquisition, J.H.L.; investigation, S.H.; methodology, S.H.; project administration, J.H.L.; software, S.H.; supervision, J.H.L.; visualization, S.H.; writing, original draft, S.H.; writing, review and editing, J.H.L.

**Funding:** This work was supported by the Advanced Biomass R&D Center (ABC) of the Global Frontier Project funded by the Ministry of Science and ICT (ABC-2011-0031354).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **An Intensified Reactive Separation Process for Bio-Jet Diesel Production**

#### **Miriam García-Sánchez 1, Mauricio Sales-Cruz 2, Teresa Lopez-Arenas 2, Tomás Viveros-García <sup>1</sup> and Eduardo S. Pérez-Cisneros 1,\***


Received: 17 June 2019; Accepted: 23 September 2019; Published: 25 September 2019

**Abstract:** An intensified three-step reaction-separation process for the production of bio-jet diesel from tryglycerides and petro-diesel mixtures is proposed. The intensified reaction-separation process considers three sequentially connected sections: (1) a triglyceride hydrolysis section with a catalytic heterogeneous reactor, which is used to convert the triglycerides of the vegetable oils into the resultant fatty acids. The separation of the pure fatty acid from glycerol and water is performed by a three-phase flash drum and two conventional distillation columns; (2) a co-hydrotreating section with a reactive distillation column used to perform simultaneously the deep hydrodesulphurisation (HDS) of petro-diesel and the hydrodeoxigenation (HDO), decarbonylation and decarboxylation of the fatty acids; and (3) an isomerization-cracking section with a hydrogenation catalytic reactor coupled with a two-phase flash drum is used to produce bio-jet diesel with the suitable fuel features required by the international standards. Intensive simulations were carried out and the effect of several operating variables of the three sections (triglyceride-water feed ratio, oleic acid-petro-diesel feed ratio, hydrogen consumption) on the global intensified process was studied and the optimal operating conditions of the intensified process for the production of bio-jet diesel were achieved.

**Keywords:** bio-jet diesel; co-hydrotreating; hydrodesulphurisation; hydrodeoxigenation; reactive distillation

#### **1. Introduction**

The area of process system engineering (PSE) has been rapidly developing since the 1950s, reflecting the remarkable growth of the oil, gas, petrochemical and biotechnological industries and their increasing economical and societal impact. It can be said that the roots of this field go back to the Sargent's pioneering school in United Kingdom [1,2]. Modelling, simulation and optimization (MSO) of large-scale (product or process) systems is a core technology to deal with the complexity and connectivity of chemical processes and their products at multiple scales [3,4]. These technologies have been implemented in easy-to-use software systems that allow the systematic solution for problem solving practitioners. The systematic (explicit or implicit) generation, evaluation and integration of a comprehensive set of design alternatives is considered a key concern for the optimal design. This systematic integration associates the PSE strongly with its traditional focus in complete plants for both, *process intensification* [5,6] and to *chemical product design* [7]. Specifically, process intensification involves the development of novel apparatuses and techniques that, in comparison with those commonly used, are expected to bring enhancements in manufacturing and processing, substantially

reducing the equipment-size/production-capacity ratio, energy consumption, or waste production, and ultimately resulting in cheaper, sustainable technologies [5]. Also, it is known that the whole scope of process intensification generally can be divided into two areas: (i) *process-intensifying equipment*, such as novel reactors, and exhaustive mixing, heat and mass transfer devices and, (ii) *process-intensifying methods*, such as new or hybrid separations, integration of reaction and separation, heat exchange or phase transition (i.e., reactive distillation), techniques using alternative energy sources (light, ultrasound, etc.), and new process-control methods (such as intentional unsteady-state operation).

The concerns about energy demand are obliging the oil-based fuels consumer countries to reconsider their energy policies to promote the investigation on trustworthy alternatives to conventional fuels. Thus, the bio-jet diesel has arisen as an alternative for petro-diesel jet fuels used in the aviation enterprises. Specifically, for the jet fuel, the International Air Transport Association (IATA) estimated that the consumption of the jet diesel would increase every year by 5% till 2030 [8]. Also, due to the growing of the flight demand and the strong regulations to diminish the CO2, IATA established a carbon neutral reduction up to 50% by 2050. In this way, the bio-jet fuel obtained from vegetable oils of from mixtures of vegetable oils and petro-diesel, can be contemplated as one of the most favourable solutions to satisfy the global demand. To now, there has been identified four main routes to obtain the bio-jet diesel: (i) oil-to-jet (deoxygenation of triglyceride and consequent hydrocracking), [9] (ii) gas-to-jet (gasification/Fischer-Tropsch reaction followed by hydrocracking), [10] (iii) alcohol-to-jet (dehydration of alcohols and successive oligomerization), [11] and (iv) sugar-to-jet (several catalytic conversions of sugars) [12]. Table 1 shows the key conversion steps, the catalyst used, the companies producing the jet-fuel and the feedstocks considered for each route to obtain bio-jet fuel.

For the conversion of oil to bio-jet fuel (OTJ), different type of oil feedstock has various converting pathways. The common pathways include the hydrogenated esters and fatty acids (HEFA) and catalytic hydrothermolysis (CH). The feedstocks for HEFA are non-edible vegetable oils, used cooking oil, and animal fats. While the feedstocks for CH are algal oils or oil plant. HEFA is a process to hydrotreat the triglycerides, saturated or unsaturated fatty acids in the non-edible vegetable oils, used cooking oils and animal fats to produce bio-jet fuels. The process is generally divided into two steps. The first step is converting unsaturated fatty acids and triglycerides into saturated fatty acid by catalytic hydrogenation, the triglycerides occur a β-hydrogen elimination reaction to yield a fatty acid during the process. The saturated fatty acid is converted to C15–C18 straight chain alkanes by hydrodeoxygenation and decarboxylation. The co-products are propane, H2O, CO and CO2. The early-developed catalysts for this step are noble metals supported with zeolites or oxides, and later shifted to other transition metals, such as Ni, Mo, Co, Mo or their supported bimetallic composites due to catalyst deactivation by poisoning, production of cracking species and process costs. The second step of HEFA is the cracking and isomerization reactions: the deoxygenated straight chain alkanes are further selectively hydrocracked and deep isomerized to generate highly branched alkanes mixed liquid fuels. The common catalysts for this step are Pt, Ni or other precious metals supported by activated carbon, Al2O3, zeolite molecular sieves. The bio-jet fuels produced by HEFA, as high energy biofuels, can be directly used in flight engine even without blending. The fuel has high thermal stability, good cold flow behaviours, high cetane number, and low tailpipe emissions, while has low aromatic content, which would cause fuel low lubricity and fuel leakage problems. Another pathway to convert algal or oil plant to jet fuel is CH, which is also named hydrothermal liquefaction (HTL). The integrated process has three steps, including pretreatment of triglycerides, CH conversion and post-refining steps. The pretreatments including conjugation, cyclization, and cross-linking, which are aiming to improve the molecular structures. The products undergo a cracking and hydrolysis reaction with the help of water and catalysts. Then it occurs catalytic decarboxylation and dehydration during the CH process. Last, post-refining hydrotreating and fractionation are designed to convert straight-chain, branched and cyclo-olefins into alkanes. The conversion process of lignocellulosic biomass to jet-fuel have advantages in lowering cost, feasible availability and no competition with food supplies. Hydroprocessed depolymerized cellulosic jet (HDCJ) is an oil upgrading technology

to convert bio-oils produced from the pyrolysis or hydrothermal of the lignocellulose into a jet fuel by hydrotreating. The main technology for bio-oil upgrading is a two-step hydroprocessing. First, the bio-oil is hydrotreated with the help of catalyst under mild conditions. Organic could be used to promote hydrodeoxygenation of bio-oil and overcome coke formation. Second, conventional hydrogenation setup and catalyst were used under high temperature for obtaining hydrocarbon fuel. The HDCJ process could produce high aromatic content, low oxygen content and few impurities jet fuel. However, there is high hydrogen consumption and deoxygenation requirements in this process, which can make a considerable expense. Moreover, the short catalyst lifetime and modest hydrocarbon yields can be challenges for used in aviation.

For the conversion of gas to bio-jet fuel (GTJ) the Fischer-Tropsch (FT) process has been commercially implemented. Fischer-Tropsch (FT) is a process to produce liquid hydrocarbon fuels from syngas. The common process for FT could be divided into six procedures: feedstock pretreatments, biomass gasification, gas conditioning, acid gas removal, FT synthesis and syncrude refining. The FT synthesis can also be divided into high temperature FT and low temperature FT. The temperature for high temperature FT is around 310–340 ◦C, and the products are main gasoline, solvent oil and olefins; the temperature for low temperature FT is around 210–260 ◦C, the products are main kerosene, diesel oil, lubricating base oil and naphtha fractions. Too low temperature of FT will format high quantities of methane as a by-product. Typical pressures of FT process are in the range of one to several tens of atmospheres. The high pressures will result the formation of long-chain alkanes. The FT fuel is free of sulphur, nitrogen, has high specific energy, high thermal stability, and cause low emissions when used for aviation. However, the disadvantages for the fuel are low aromatic content and less energy density, which would also cause a low efficiency and high production cost for the process.

For the conversion of sugar to bio-jet fuel (STJ) the biotechnological process of convert sugars to alkane-type fuels directly instead of firstly converting to ethanol intermediate, which is called Direct Sugar to Hydrocarbons (DSHC) has been implemented commercially (see Table 1). The feedstock for the DSHC are similar to the feedstock of ethanol production, including the sugar cane, beets and maize. DSHC is a process to produces alkane-type fuels directly from sugars via fermentation. It is different from the alcohol to jet pathway, which needs an alcohol intermediate. The technology is developed based on the development of genetic engineering and screening technologies that enable to modify the way microbes metabolize sugar. A complete conversion process of DSHC, involving six major steps: pretreatment and conditioning, enzymatic hydrolysis, hydrolysate clarification, biological conversion, product purification and hydroprocessing has been proposed [13]. The DSHC has a low energy input due to the low temperature of the fermentation, while the fuel blend is limited (10%) and not meet some performances standards and it is also identified as more suitable for production high-value chemicals.

Finally, the conversion route to convert alcohols to bio-jet fuel (ATJ) use several processes depending of the feedstock. The process of production hydrocarbons in the jet fuel range from the alcohols generally undergoes a four-step upgrading process. First is the alcohol dehydration to generate olefins, then the olefins are oligomerised in the presence of catalysts to produce a middle distillate. Next, the middle distillates are hydrogenated to produce the jet-fuel-ranged hydrocarbons and a final step is the distillation to purify the bio-jet fuel product. Commercial production always use ethanol, butanol and isobutanol to be the intermediate to converse biomass to jet fuel. The economics of alcohol to jet process is mainly affected by the way to produce alcohol, while the biochemistry way has relatively smaller minimum jet fuel selling price (MJSP), the sugar cane and starch are suitable feedstock from an economics prospective. Complete reviews of the technologies described above considering the techno-economic analysis are given in [14–17].


**Table 1.** Summary of bio-jet fuel production pathways. Data obtained from [14].

Considering the above technologies, the (i) oil-to-jet and (ii) gas-to-jet are contemplated as the most convincing alternatives in the short term. In fact, the bio-jet fuels produced with these technologies are now allowed by ASTM specification D7566-14 for blending into commercial jet fuel at concentrations up to 50% [18]. Recently [19], a hydrotreating reactive distillation column (RDC) for the production of green-diesel from sulphured petro-diesel and non-edible vegetable oils mixtures was proposed. From this work, it was concluded that the deep hydrodesulphurisation (HDS) of petro-diesel is strongly affected by feeding high amounts of triglycerides to the RDC, while it is not affected when only fatty acids are fed. Also, an integrated reactive separation process for the production of jet-diesel was developed but the assumptions about the conversion of tryglycerides to fatty acids and the yield of isomerization-cracking of the linear hydrocarbon chains lead to an oversimplified reactive separation process [20]. Therefore, in this work, from a process intensification context, an innovative intensified three-step reactive separation process for bio-jet diesel production is proposed.

Thus, the intensified three-step reactive separation process consist of: (1) a triglyceride hydrolysis section where a catalytic heterogeneous reactor is used to convert the triglycerides to the resultant fatty acids, followed by a three-phase separation device and a sequence of two distillation columns to obtain pure fatty acid and glycerol; (2) a hydrotreating section with a reactive distillation column used to simultaneously carried out the deep hydrodesulphurisation (HDS) of petro-diesel and the hydrodeoxigenation, decarbonylation and decarboxylation of the fatty acids; and (3) an isomerization-cracking section with a hydrogenation fixed bed reactor connected to a two phase flash separator to produce bio-jet diesel with the required fuel properties. It should be pointed out that the conversion and yield assumptions used in [20] are not considered here and the appropriate reaction kinetics found in the open literature [21,22] is used to perform the rigorous simulation of the intensified process.

#### **2. The Intensified Reactive Separation Process**

The production of bio-jet fuel can be described by the subsequent reactions (see Figure 1): (i) hydrogenation of the C=C bonds present in the tryglycerides, (ii) hydrogenolysis of the tryglycerides to

produce the respective fatty acids, (iii) deoxygenation of the resultant fatty acids to obtain n-paraffins and (iv) hydroisomerization and hydrocracking of the n-paraffins to produce a mixture containing mainly isomerized shorter chains (*i*-C8−*i*-C16) that are appropriate as jet fuel. The first three reactions can be carried out using a metal catalyst or a MoS2-based catalyst. For the hydroisomerization and hydrocracking reaction (iv), a metal/acid bifunctional catalyst is recommended. Figure 1 shows the reaction pathways for the hydroconversion of tryglycerides into bio-jet fuel. It has been shown [23] that the direct hydrogenation of tryglycerides may be very expensive due to the price of the hydrogen at high pressures. Thus, the conventional hydrolysis of tryglycerides with high pressure steam is preferred. An alternative to the hydrogenation and hydrolysis at high pressure is to use a heterogeneous hydrolysis catalytic reactor [24] at moderated pressures. Thus, in the present work this technological alternative is considered.

**Figure 1.** Reaction pathways of the hydroconversion of tryglycerides into biojet fuel.

Figure 2 shows a simplified flow sheet of the intensified three-step reactive separation process for bio-jet diesel production. The first section (hydrolysis) consists of two pumps operating from 1 to 30 atm and two heat exchangers where the exit temperature is set to 280 ◦C. The heterogeneous hydrolysis reactor operates at 30 atm and it consists of two 25 m length tubes with 0.5 m of diameter. Further, a sequence of a three-phase flash and a two-phase drum are used to separate the fatty acid-glycerol-water mixture and the unconverted glycerides that are recycled to the hydrolysis reactor. Finally, two distillation columns are used to separate the pure fatty acid from a low concentration glycerol-water mixture and to produce pure glycerol. The second section, hydrodesulphurisation-hydrodeoxigenation (HDS-HDO) consists of a reactive distillation column (RDC) where a mixture of petro-diesel and fatty acid is fed to the RD column at stage 9 and excess hydrogen is fed at the bottom. The gases produced by the HDS-HDO reactions are released from the two-phase condenser. The exit liquid mixture of C11–C12 linear hydrocarbon chains is mixed with the exit stream containing C17–C18 hydrocarbons produced by the HDO reactions and the ultra-clean petro-diesel products. The third section consists of a heterogeneous hydro-conversion reactor where the long linear hydrocarbon chains (C11–C18) are isomerized and cracked to attain the bio-jet diesel.

**Figure 2.** Simplified three-step reaction-separation process for bio-jet diesel production.

#### *2.1. The Triglycerides Hydrolysis Section*

Vegetable oils and fats have been considered as one of the most used renewable raw materials in the chemical industry. These can be hydrolyzed to produce free fatty acids (FFA) with a high degree of purity to be used in the synthesis of chemically pure compounds [25]. However, in the present work, the proposed technological alternative to produce bio-jet fuel considers non-edible oils, waste oils and animal fats as feedstock, which are not in conflict with food resources. Fatty acids are used in a wide variety of industries, for example in the pharmaceutical and cosmetics industry. Also, fatty acids can be utilised to produce n-alkane chains through a decarboxylation process [26]. These hydrocarbons work properly in internal combustion engines as substitutes for petro-diesel. Actually, the hydrolysis of fats and vegetable oils, composed mainly of triglycerides, has been practiced in the industry for many years. In general, the hydrolysis of the esters occurs through the acyl-oxygen break [27] with an excess of water at high temperature or using an appropriate acid catalyst to hydrolyze the glycerol pillar in the ester group of any triglyceride (TG), diglyceride (DG) or monoglyceride (MG) [28]. The result of the hydrolisis reactions is the production of three moles of FFA and one mole of glycerol (Gly). In this work, triolein is considered as the main triglyceride compound and oleic acid as the correspondent fatty acid. The three consecutive reversible reactions can be written as:

$$\text{triolein} + \text{H}\_2\text{O} \Leftrightarrow \text{diolein} + \text{oleic acid} \tag{1}$$

$$\text{iodolein} + \text{H}\_2\text{O} \Leftrightarrow \text{mono-olein} + \text{oleic acid} \tag{2}$$

$$\text{monno-olein} + \text{H}\_2\text{O} \Leftrightarrow \text{glycerol} + \text{oleic acid} \tag{3}$$

Commonly, the pure fatty acids are obtained from the reaction of vegetable oils and/or animal fats with superheated steam. The commercial hydrolisis conditions are around 100–260 ◦C and 10–7000 kPa using a 0.4–1.5 wt % water-to-oil ratio. Several variants of this technology have been used by industry [29,30]. In this work, the hydrolisis kinetics reported in the open literature [21] is used for the numerical simulation of the catalytic reactor using tungstated zirconia (WZ) and the solid composite SAC-13 as catalyst.

#### *2.2. The HDS-HDO Reactive Distillation Section*

The hydrotreatment of non-edible vegetable oils, waste oils or animal fats, this is, oils and fats that are not used for food and others medical applications, to produce renewable fuels has several advantages: (i) flexibility in the disposal of raw materials due to the great variety of oils available from vegetables on the earth; (ii) the process can be carried out using the existing infrastructure in petro-refineries and (iii) the bio-fuels produced can be used in conventional internal combustion machines since these bio-fuels have properties similar to those obtained from mineral oil [31]. The hydrotreatment of vegetable oils can be accomplished using traditional catalysts, for example with NiMo/Al2O3 catalysts. The hydrotreatment mainly produces n-paraffins through the hydrodecarboxylation, hydrodecarbonylation and hydrodeoxigenation reactions. From the refining point of view, the hydrodecarboxylation reaction is better than the hydrodeoxygenation reaction since it consumes less hydrogen. However, the large amount of CO (or CO2) generated represents a problem for refining, since CO2 can form carbonic acid with liquid water [32]. This means that the risk of carbonic corrosion of the reactive separation equipment is a key design problem. In general, vegetable oils can be hydrotreated as pure compounds or can be mixed with petro-diesel to be co-hydrotreated, in such a way that the hydrodesulphurisation (HDS) and hydrodeoxygenation (HDO) reactions are carried out simultaneously in a single unit. Therefore, the reactions considered for the HDS of the sulphured petro-diesel can be written in a simplified form as:

$$\text{Thiopphene (Th)} + 2\text{ H}\_2 \rightarrow \text{Butadiene} + \text{H}\_2\text{S} \tag{4}$$

$$\text{Benzothi} \\ \text{sphere (BT)} + 3 \text{ H}\_2 \rightarrow \text{Ethylbenzene} + \text{H}\_2\text{S} \tag{5}$$

$$\text{DBT} + 2\text{ H}\_2 \rightarrow \text{Biphenyl} + \text{H}\_2\text{S} \tag{6}$$

For the HDS reactions only the hydrogenolysis reaction pathway is considered (Equation (6)). The simplified HDO reactions of the oleic acid can be written as:

$$\text{R}\text{ }\text{C}\_{18}\text{H}\_{34}\text{O}\_{2} + 4\text{ H}\_{2} \rightarrow \text{n-C}\_{18}\text{H}\_{38} + 2\text{ H}\_{2}\text{O} \quad \text{(hydrodeoxygmentation)}\tag{7}$$

$$\text{C}\_{18}\text{H}\_{34}\text{O}\_2 + 2\text{ H}\_2 \to n\text{-C}\_{17}\text{H}\_{36} + \text{H}\_2\text{O} + \text{CO} \quad \text{(decarbonym)}\tag{8}$$

$$\rm C\_{18}H\_{34}O\_2 + H\_2 \to n-C\_{17}H\_{36} + CO\_2 \quad \text{(decarbohydation)}\tag{9}$$

The reaction rate equations for the HDS reactions can be obtained from [19] and, due to the absence of reliable kinetic expressions for the complex HDO reactions, a 99% conversion of oleic acid is assumed.

#### *2.3. The Isomerization-Hydrocracking Section*

In general the processes of isomerization and cracking through hydrogenation occur simultaneously. When hydrogenation favours isomerization it is called hydro-isomerization and when it promotes cracking it is called hydro-cracking. Depending on the characteristics of the hydrocarbon chains to be hydrotreated, the selection of the appropriate catalyst for this purpose can be made in a wide range of possibilities [33]. Specifically, for the production of jet diesel, the catalysts are bifunctional and are characterized by having acid sites that allow the selective function of isomerization and cracking simultaneously. Experimental evidence [34] indicate that the conversion of the n-paraffins can be described by the following reactions network:

where each single reaction is irreversible and follows a pseudo first-order kinetic. (A) is the concentration of the n-paraffin; (B) is the concentration of the iso-paraffin; (C), (D) and (E) are the concentrations of the cracking products. Following the linear reaction sequence path: A–>B–>C the hydro-conversion reactions could be written as:

$$n\text{-C}\_{11}\text{H}\_{24} \rightarrow i\text{-C}\_{11}\text{H}\_{24} \tag{10}$$

$$\text{m-C}\_{12}\text{H}\_{26} \rightarrow \text{i-C}\_{12}\text{H}\_{26} \tag{11}$$

$$\rm{m-C}\_{13}\rm{H}\_{28} \rightarrow \rm{i-C}\_{13}\rm{H}\_{28} \tag{12}$$

$$n\text{-C}\_{14}\text{H}\_{30} \rightarrow i\text{-C}\_{14}\text{H}\_{30} \tag{13}$$

$$n\text{-C}\_{16}\text{H}\_{34} \rightarrow i\text{-C}\_{16}\text{H}\_{34} \tag{14}$$

$$\text{i-C}\_{17}\text{H}\_{36} \rightarrow \text{i-C}\_{17}\text{H}\_{36} + \text{H}\_2 \rightarrow \text{i-C}\_8\text{H}\_{18} + \text{i-C}\_9\text{H}\_{20} \tag{15}$$

$$\rm{i-n-C\_{18}H\_{38}} \rightarrow \rm{i-C\_{18}H\_{38}} + \rm{H\_2} \rightarrow \rm{i-C\_8H\_{18}} + \rm{i-C\_{10}H\_{22}} \tag{16}$$

The kinetic equations for the isomerization and cracking of the linear hydrocarbon chains were taken from [16] and a catalyst of Pt supported on a nano-crystalline large-pore BEA zeolite is used [28].

#### **3. Results and Discussion**

In order to determine the best operating conditions of the intensified three-step hydrotreating reactive-separation process, the effect of different design and operating variables on the performance of the global process was investigated and analysed through intensive simulations. The intensive simulations of the intensified hydrotreating reactive separation process were carried out in the commercial Aspen-Plus simulator environment. Thus, heat and mass transfer phenomena and mixing issues were not taken into account.

#### *3.1. Hydrolysis Section*

Different amounts of triolein-water feed ratio were used for the simulation of the intensified process. The numerical results shown in Table 2 are for the reference case (production of 70 kmol/h of pure oleic acid). From such numerical results it is found that total conversion of triolein is attained at 553 K and excess water (265 kmol/h). This water flow corresponds to 1/9 triolein/water feed ratio. The three-phase flash that separates water, oleic acid-water-glycerol and triolein-diolein mixture operates at 5 atm and 410 K, and the flash drum for the separation of the water-oleic acid-glycerol-unconverted glycerides operates at 1 atm. The two distillation columns to produce pure oleic acid and pure glycerol operate at 1 atm. It should be pointed out that the non-ideality of the reactive polar mixture is modelled using the RK-ASPEN equation of state.


**Table 2.** Simulation results for the reference case (hydrolysis section).

#### Effect of the Triolein-Water Feed Ratio

The effect of the triolein-water feed ratio was studied in order to verify the maximum conversion of triolein to oleic acid. Figure 3a shows the effect of the feed ratio on the hydrolysis reactor exit composition. It can be observed that as the feed ratio increases, the triolein composition decreases to zero at 1/7 triolein-water feed ratio, approximately. However, it should be noted that, as the feed ratio increases, the hydrolysis reactor exit mixture contains less fatty acid (oleic acid). This is expected due to the excess of water fed to the reactor. Also, it can be observed that after a feed ratio of 1/7, the final reactor exit mixture contains only water, oleic acid and glycerol. Figure 3b shows the effect of the triolein-water feed ratio on the distillation column bottom flow (stream 11 in Figure 2). It can be noted in Figure 3b that as the feed ratio increases, the amount of oleic acid produced increases to a constant value of 88 (kmol/h). The production of 70 (kmol/h) at the bottom of the distillation column at a feed ratio of 1/9 has been taken as a reference study case.

**Figure 3.** (**a**) Effect of the triolein-water feed ratio on the hydrolysis reactor exit composition. (**b**) Effect of the triolein-water feed ratio on the distillation column bottom flow (stream 11).

#### *3.2. HDS-HDO Section*

In order to perform the deep hydrodesulphurisation of the sulphured petro-diesel and the reactions involved in the hydrodeoxigenation of the fatty acid, a two-zone reactive distillation column (RDC) is recommended [13]. Figure 4 shows the sizing details of the RDC, as well as the basic information for the simulation of the reference case. The RDC consists of 14 stages with two sections of reactive stages (5–7, 9–11) operating at 30 atm. It can be observed in Figure 4 that the liquid streams C11–C12 and C13–C18 contains mainly the larger linear hydrocarbon chains. Table 3 displays the numerical simulation results for the HDS-HDO section with 70 kmol/h of oleic acid and 100 (kmol/h) of petro-diesel fed at stage 9. Hydrogen feed was set to 400 (kmol/h) in order to accomplish full conversion of oleic acid and deep HDS of petro-diesel.


**Table 3.** Simulation results of the HDS-HDO section in ASPEN-PLUS (reference case).

**Figure 4.** Hydrotreating two-zone reactive distillation column. Oleic acid premixed with sulphured petro-diesel fed at stage 9.

#### 3.2.1. Effect of the Oleic Acid-Petro-Diesel Feed Ratio

The temperature profile along the reactive distillation column (RDC) assuming that oleic acid is premixed with the petro-diesel and fed in stage 9 is shown in Figure 5a. The co-hydrotreating process starts with a low concentration of oleic acid (10 kmol/h) in the mixture fed and it was continuously augmented up to 70 (kmol/h). It can be seen in Figure 5a that as the content of oleic acid increases, the temperature at the reactive zone I is reduced. The temperature reduction can be explained by considering that the boiling point of oleic acid is higher than the boiling point of most hydrocarbons present in the sulphured petro-diesel. The effect of the oleic acid composition in the feed on the deep HDS of the sulphured petro-diesel is shown in Figure 5b. It can be noted in Figure 5b that, even for a high content of oleic acid (70 kmol/h), the deep HDS of petro-diesel is accomplished. The analysis of these results leads to the conclusion that the performance of the reactive distillation column (RDC), initially used for deep hydrodesulphurisation of petro-diesel, is not affected by the introduction of a vegetable oil with a high content of fatty acid (oleic acid).

**Figure 5.** (**a**) Temperature profile along the RDC. (**b**) DBT liquid composition profile.

#### 3.2.2. Hydrogen Consumption and Liquid Water Production

The molar flow of hydrogen required for the HDS and HDO reactions as a function of the oleic acid content in the feed is shown in Figure 6a. It can be noted in Figure 6a that, for an oleic acid feed flow greater than 10 (kmol/h), the temperature of the condenser is reduced to 215 ◦C, and for a feed flow of 70 (kmol/h), the required hydrogen is increased to 400 (kmol/h) and the temperature of the condenser should be diminished to 185 ◦C. It should be pointed out that, if the temperature of the condenser is not decreased, the numerical convergence of the RDC simulation is not reached. Thus, a control loop between the source of vegetable oil (oleic acid) and the hydrogen supplied should be linked to the temperature of the condenser. The liquid water composition profiles along the RDC for different oleic acid feed flows are shown in Figure 6b. It can be noted in Figure 6b that the concentration of liquid water increases from reactive zone I (steps 5–7) to the upper part of the RDC. Liquid water is mainly produced in reactive zone II (reactive stage 9–11) by the hydrodeoxygenation and decarbonylation reactions and its molar fraction increases continuously from stage 9 to the upper part of the RDC, as the concentration of oleic acid in the feed is augmented. Therefore, in order to minimize catalyst deactivation and corrosion of the RDC equipment, the amount of oleic acid in the mixture fed should not be higher than 70 (kmol/h).

**Figure 6.** (**a**) Hydrogen feed flow required on stage 12. Oleic acid-petro-diesel mixture fed at stage 9. (**b**) Water liquid composition profile.

#### 3.2.3. Hydrocarbon Distribution and Release of Generated Gases

Figure 7a shows the hydrocarbon liquid composition distribution. It can be noted in Figure 7a that the bottom exit stream contains mainly the linear hydrocarbons C17, C18, C16, C14, and C13. At the top of the RDC a rich liquid mixture of C11 and C12 linear hydrocarbons is obtained. It should be pointed out that such linear hydrocarbons (C11–C18) are further mixed and fed to the isomerization-hydrocracking reactor (see Figure 2). Figure 7b shows the gas composition profile along the RDC. In the case of the hydrotreating RDC, the gases produced are mainly driven in the vapour phase at each equilibrium stage and released from the partial condenser at the top of the RDC. It can be noted in Figure 7b that the gases produced by HDS and HDO reactions: CO2, CO, steam and H2S and the unconverted excess of hydrogen are released at the top of the column. It is interesting to observe that the vapour composition of C11 and C12 decreases at the top of the column.

**Figure 7.** (**a**) RDC profiles for oleic acid-petro-diesel blend fed at stage 9 at 30 atm. (**b**) Gas compositions profile along the RDC. Reproduced with permission from García-Sánchez et al., Proceedings of the 28th European Symposium on Computer Aided Process Engineering; Published by Elsevier B.V., 2018. [20].

#### 3.2.4. Operability and Controllability of the HDO-HDS Reactive Distillation Section

Reactive distillation column mathematical models are highly nonlinear, and multiple steady states (MSS) solutions have been reported by many researchers [35–37]. However, none of these works have addressed the HDO-HDS reactive distillation process, most of them studied the MTBE and TAME cases. In the present work, the MSS solutions are only briefly described. Two multiplicity types can be found: input and output multiplicities, but a combined input–output multiplicity may also be present. Output multiplicity occurs when one set of input variables (manipulated variables) results in two or more unique and independent sets of output variables (measured variables). In chemical reactors and reactive distillation columns there are usually three steady states associated to ignition (high conversion), extinction (low conversion) and medium unstable conversion. On the other hand, input multiplicity occurs when two or more unique sets of input variables produce the same output condition. This type of multiplicity has important implications for close loop control since it is related to the so-called zero dynamics of the system, which is associated with unusual, unexpected, or inverse responses of the outputs after a step change, has been applied to the inputs. As RDC share some common features with chemical reactors and conventional (non-reactive) distillation columns, the RDC behavior may exhibit input–output multiplicity with three steady states or output multiplicity, with a large number of different steady states induced by ignition and extinction of every single reactive column tray [38]. Thus, in the present work, from the intensive simulation results, it was found that the reflux ratio is a key parameter for an appropriate RDC design and operation. Additionally, an increment in the hydrocarbon feed (sulphured diesel) flowrate leads to lower conversions of the organo-sulfur and fatty acid compounds with the reduction of the H2/HC feed ratio. In addition, it was found that, the DBT conversion was highly affected by the variations of the reboiler heat duty, while the fatty acid conversion was practically constant (~99%). However, DBT conversion design target (99.9%) was achieved at four different reboiler duties, indicating the existence of input-output multiplicities. It was also determined that the amount of catalyst loaded at the reactive zones must be greater than 8000 kg in order to achieve an ultra low sulfphur diesel (ULSD) production and a complete conversion of the fatty acids. This is, if a suitable excess of H2 is present in the reaction zones, the behavior of the RDC is very likely to a HDS conventional catalytic reactor. Therefore, it may be said that the multiplicities found in the intensive simulations are highly related to the specific phenomena (reaction-separation) involved. It should be pointed out that the accurate determination of the thermodynamic properties of all species participating in the reacting mixture must be taken carefully. The correct determination of the boiling points, critical properties, etc. of all species involved in the reactive separation process and

the computation of the complex phase equilibrium in the intensified reactive separation process is a key step. This is, for example, for the HDO-HDS section, the inaccuracy of the property values (boiling point of the fatty acids) can lead to multiple steady states and unstable operation and control of the reactive distillation process. The analysis of the intensive simulation results suggest that, under optimal design and operating conditions, reactive distillation can be considered as a viable technological alternative to produce bio-jet diesel trough a RDC co-hydrotreating process.

#### *3.3. The Isomerization-Cracking Section*

The numerical simulation results with the final composition of the bio-jet diesel obtained in the hydro-conversion reactor are shown in Table 4. It should be mentioned that an excess of hydrogen must be fed to the reactor in turn to accomplish the complete isomerization and cracking of the longer hydrocarbon chains. The total pressure of the reactor is set at 80 atm and a 5:1 hydrogen-hydrocarbon feed ratio was considered. The exit stream from the isomerization-cracking reactor is further flashed to eliminate the hydrogen excess and the light hydrocarbons produced during the cracking step. The isomerization-cracking kinetics and reactor arrangement were taken from the open literature [34]. From Table 3 it can be observed that the bio-jet-diesel produced contains mainly *i*-C8 to *i*-C16 hydrocarbons indicating that a light bio-jet diesel is produced.


**Table 4.** Simulation results of the isomerization-cracking section in ASPEN-PLUS (reference case).

#### **4. Conclusions**

An intensified three-step hydrotreating reaction-separation process for the production of bio-jet diesel from triolein and petro-diesel mixtures has been developed. Through intensive simulations the effect of different operating variables (triglyceride-water feed ratio, oleic acid-petro-diesel feed ratio, hydrogen consumption) on the performance of the intensified reactive separation process was studied. By analysing the simulation results, it can be established that, a 1/9 triolein-water feed ratio guarantee the complete conversion of triolein to oleic acid at moderated pressures (30 atm) at the hydrolysis section. For the HDS-HDO section, it can be mentioned that a mixture containing up to 70 (kmol/h) of oleic acid in the hydrocarbon feed of the RDC is the optimal mixture composition in order to achieve the deep hydrodesulphurisation of petro-diesel. Also, with this mixture composition, any change to the basic structure (reactive and non-reactive stages) of the RDC is not required. It should be pointed out that the HDS-HDO RD column should operate at moderated pressures (30 atm) in

order to allow the fast vaporization of the light gases produced and avoid the corrosion problem with carbonic acid generated by the reaction of CO or CO2 with water. For the isomerization and cracking section, a low-pressure flash separator, after the hydro-conversion reactor is required to eliminate the undesirable side products. Therefore, it can be concluded that the key design and operating parameters for the production of the bio-jet diesel are: (i) in the hydrolysis section, the water excess and the total pressure for the heterogeneous catalytic hydrolysis reactor; (ii) for the HDS-HDO section, if high molar flows of fatty acid are considered, it is mandatory to have more reactive stages in the HDS-HDO reactive distillation column in order to achieve ultra-clean (no-sulphur) petro-diesel at the bottom of the column. Finally, the absence of light compounds in the exit stream of the heterogeneous isomerization-cracking reactor is required in order to circumvent undesirable side products and reach the suitable fuel properties required by the international standards. Finally, the economic and sustainability analysis of the intensified reactive separation process to produce bio-jet diesel is being carried out.

**Author Contributions:** Conceptualization, M.S.-C. and E.S.P.-C.; Formal analysis, T.V.-G.; Investigation, E.S.P.-C.; Methodology, M.G.-S. and E.S.P.-C.; Software, T.L.-A.; Validation, M.S.-C. and T.V.-G.; Writing—original draft, M.G.-S.; Writing—review & editing, T.L.-A. and E.S.P.-C.

**Funding:** This research was funded by CONACyT PhD scholarship of Miriam García-Sánchez.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Journey Making: Applying PSE Principles to Complex Curriculum Designs**

#### **Ian Cameron \* and Greg Birkett**

School of Chemical Engineering, The University of Queensland, Brisbane 4072, Australia; g.birkett@uq.edu.au **\*** Correspondence: itc@uq.edu.au

Received: 7 November 2019; Accepted: 16 March 2020; Published: 23 March 2020

**Abstract:** Since the 1950s, Process Systems Engineering (PSE) concepts have traditionally been applied to the process industries, with great effect and with significant benefit. However, the same general approaches and principles in designing complex process designs can be applied to the design of higher education (HE) curricula. Curricula represent intended learning journeys, these being similar to the design of process flowsheets. In this paper, we set out the formal framework and concepts that underlie the challenges in design of curricula. The approaches use generic and fundamental concepts that can be applied by any discipline to curriculum design. We show how integration of discipline-specific concepts, across time and space, can be combined through design choices, to create learning journeys for students. These concepts are captured within a web-based design tool that permits wide choices for designers to build innovative curricula. The importance of visualization of curricula is discussed and illustrated, using a range of tools that permit insight into the nature of the designs. The framework and tool presented in this paper have been widely used across many disciplines, such as science, engineering, nursing, philosophy and pharmacy. As a special issue in memory of Professor Roger W.H. Sargent; we show these new developments in curriculum design are similar to the development of process flowsheets. Professor Sargent was not only an eminent research leader and pioneer, but an influential educator who gave rise to a new area in Chemical Engineering, influencing its many directions for more than 50 years.

**Keywords:** process systems engineering; design; higher education; curricula; visualization

#### **1. Introduction**

Curriculum design stands at the heart of all education. It is a multiscale challenge across different time and space scales—from whole-of-curriculum design considerations, to distinct learning units or modules, down to the day-to-day learning elements at the lowest level of consideration. It spans sequential stages of learning—from early learners, primary, secondary, tertiary and continuing professional development (CPD). By necessity, curriculum designs seek to embody stated intended outcomes for learners that address knowledge domains, application of knowledge, and personal- and professional-attribute development.

It might seem strange to some that curriculum design could be intimately related to Process Systems Engineering (PSE) thinking and application. In what follows, we show the development of complex curricula from the basic underlying concepts and building blocks that mirror many aspects of PSE. In doing so, we emphasize that PSE possesses a much broader interpretation and application than has traditionally been adopted.

In 1967, Roger Sargent wrote in *Chemical Engineering Progress* a review on "Integrated Design and Optimization of Processes". He stated the following: "Although we are in sight of a truly integrated approach to design of complete processes, a great deal of work remains to be done. With the need for more sophisticated analysis of larger complexes, it is more important than ever to join hands with those working in the fields of control engineering, operational research, numerical analysis and computer science" [1].

The discipline area of PSE arose from the application of systems engineering concepts to industrial processes [2,3]. The 1960s was a period of rapid digitalization in industry, affording significant advances in modeling, control, optimization and new computer-based numerical methods. The focus on "engineering" of "systems" that were primarily within the "process" industries drove PSE as a new focus within Chemical Engineering. The concept of "engineering" as ingenuity in design, using thinking and practices around a set of things that work together ('system'), can be applied to the general idea of any "process". That liberates PSE from the narrow confines of industrial and manufacturing sectors.

The 1967 statements of Roger Sargent, adapted to the case of curricular design, ring true. The tasks require an integrated approach that ensures the final curriculum design is "fit-for-purpose". It is a complex set of tasks dealing with many interconnected learning units, their attributes and intended outcomes. The appropriate sequencing of learning, as well as generating deep insights into the nature and behavior of the design, is essential. It is also a task where the skills and insights of numerous people are necessary to arrive at designs that deliver the requisite outcomes. In short, it bears many resemblances to traditional PSE thinking and practices.

Curriculum design practice has a long and important history. The rapid expansion of human knowledge in all professional domains has increased the need for learning designs that must meet the demand of current and future work demands in ever-changing environments. That is a long-standing challenge, and one that continues to challenge educators.

Early work by Dewey [4] and others, such as Tyler [5], set the scene for modern curriculum considerations. The expanding digitalization trends across society, with the creation and growth of the internet, have brought the need to use enabling information technologies, visualizations and user-centered web systems, to improve curriculum designs and their deployment. These information and communication technologies (ICT) can enhance the curriculum design process. The design process relies on structured information and its use to create educational pathways for learners.

In what follows, we outline why curriculum design and deployment is important, the basic concepts and processes that help engineer the learning system and our ability to assess designs in both qualitative and quantitative ways. Section 2 deals with design purpose and practice, before turning in Section 3 to the fundamental concepts and building blocks that permit designers to assemble a range of desired learning pathways. Section 4 shows the principles in practice, using an actual case study in Chemical Engineering at The University of Queensland. The final section considers how many other disciplines have derived benefit from a systems design environment that was pioneered within Chemical Engineering.

#### **2. Curriculum Design: Purpose and Practice**

#### *2.1. Purpose*

Curriculum designs set out the key learning pathways that an educational organization places before a learner. That is, curricula represent learning journeys, hence the title of this paper: "Journey Making".

"Curriculum" is related to the Latin word "currere", meaning "to run". As such, it speaks of a *pathway* or *a course* traversed by participants to reach a goal. End-goals and intermediate goals are crucial outcomes within curricula. The attainment of such goals is incremental.

In the last few decades, the underlying principles of curricula have been revisited and emphasized by prominent educational researchers across many disciplines [6,7]. These researchers and practitioners have enhanced earlier understandings of the role of curricula, the educational psychology around learning and the science of learning [8,9].

Barnett and Coate [6] recognized "that curricula have distinctive, but integrative components, as well as allowing for different weightings of each domain within any one curriculum" (p. 70).

It is not just academic researchers who are interested in curriculum design and outcomes. Due to the professional nature of engineering registration, practice and graduate mobility, global organizations such as the International Engineering Alliance (IEA), as represented by the Washington Accord [10], "establish and enforce internationally bench-marked standards for engineering education and expected competence for engineering practice".

In the case of the IEA, 29 countries that span the globe are signatories to that agreement. The professional engineering bodies or accreditation agencies in each country administer the accreditation of programs and curricula for higher education institutions that produce graduate engineers. In the USA, the Accreditation Board of Education & Technology (ABET, Baltimore, MD, USA) administers undergraduate degree programs [11]. Other organizations such as the European Network for Accreditation of Engineering Education (ENAEE), administer the EUR-ACE accreditation system for engineering graduates [12].

This means that, for accreditation of engineering programs across most of the globe, there are necessary learning outcomes that must be seen within the curriculum design, and importantly the evidence gathered to show that graduate engineers display those outcomes to an acceptable level. As well as the professional incentives for curriculum design, there are important institutional incentives in terms of developing curricula that sets apart one institution from another. This is often seen in terms of the learning pathways that students travel through during their undergraduate and graduate programs.

This raises the importance of excellence in curricular design and delivery. PSE-type approaches can provide the rigor necessary to achieve innovative, accredited curricula that provide flexibility in learning pathways to reach intended outcomes. PSE principles related to output requirements, integrated system and quality control through measurement and assessment strategies naturally fit into curriculum considerations.

#### *2.2. Practice*

Curriculum design practice across the higher education (HE) sector is extremely varied in the processes adopted to imagine what needs to be done, as well as the tools that might aid in developing, understanding and displaying designs. Many institutions struggle to properly design and document outcomes. Many designs are done by a small group of discipline experts with little 'buy-in' from colleagues. Many academics are not aware of the design principles and how their specific learning units integrate into the whole curriculum.

In recent years, a number of tools have been developed to help disciplines address the design activities in a structured manner. These primarily are mapping tools for learning outcomes (LOs), and they are limited in scope for doing comprehensive curriculum design across all higher education disciplines [13–15]. Some systems such as SOFIA have more pathway features [16]. Early comprehensive work in this space was done by the authors within Chemical Engineering [17]. PSE approaches can help in organizing curriculum around key elements that are combined to address intended learning outcomes and tracking those outcomes through the curriculum. This can help with the design process.

#### **3. Curriculum Design: Underlying Concepts, Building Blocks and Pathways**

In this section, we investigate the fundamental concepts upon which curriculum design rests. We consider the following questions:


#### *3.1. Basic Concepts*

What are we trying to achieve in curriculum design? This is a crucial question to consider before undertaking any design activities. It mirrors the focus of Requirements Engineering in seeking to define a functional specification of the needs and desires of stakeholders [18]. The general equivalent in curriculum is the Intended Learning Outcomes (ILOs). A survey of the key literature, including national and international accrediting agencies and research monographs, shows that educational program requirements are in three major areas [11,19–21]. Figure 1 displays the schema of Barnett and Coate [6] that covers these requirements.


There is an intentional overlap in this schema, emphasising that all three elements combine synergistically toward the development of graduate capabilities. This schema is not confined to engineering disciplines, but has universal application.

**Figure 1.** Educational schema as the major focus of educational designs.

Figure 2 shows the engineering program outcomes schema for Engineers Australia, the peak professional body in Australia. The accreditation frameworks in the USA and Europe, including The Washington Accord signatories, have a similar schema.

**Figure 2.** Engineers Australia, Stage 1 competency schema.

From this high-level schema, it is clearly necessary to develop more detailed outcomes at the level of the individual learning-unit level and show how integration across learning units addresses the outcomes within the *Knowing*, *Acting* and *Being* schema.

#### *3.2. Building Blocks*

Similar to flowsheet development that deals with processing units and the connected streams, curriculum development requires description of the learning units or modules. In some institutions, these are called "courses" or "subjects". They represent the lowest-level learning and teaching element in a degree program. The connections between units show required prior learnings (RPL). Besides these topological elements, the development of Intended Learning Outcomes (ILOs) for individual learning units, and subsequently the whole curriculum, require a wide range of other important elements. Figure 3 shows the key elements related to the process leading to the "curriculum as designed" stage.

**Figure 3.** Curriculum design building blocks and information.

It should be noted that, beyond "curriculum as designed", there are two further concepts: "curriculum as delivered" and "curriculum as experienced". These concepts relate to the actual delivery of the curriculum by instructors and staff to student engineers, and the curriculum as actually experienced by the students. These are important concepts, but discussion of these is beyond the focus of this paper.

It should be noted that, in many higher education institutions, students have the ability to design their own personalized curriculum. Currently, for many students, it is difficult for them to see the connections between learning unit choices, sometimes leading to fragmentation in the overall curriculum. Moreover, over time, well-designed curricula can fragment, as learning units change with little consideration given to the implication of those changes on other integrated learning units.

#### 3.2.1. Information Management

In order to gain maximum benefit from the development of a computer-aided environment for curriculum design, it is important that the information contained in these elements is both organized and extensible within data structures. We can classify the elements into two main categories:


Generic concepts are not discipline-specific, but need to be extensible and editable. They could include learning activities, attainment levels or assessment strategies. They are vital in developing the learning-outcome statements. Specific concepts relate to the discipline area considered in the design. Obvious specific areas would be discipline knowledge and application.

Taxonomies can be used to capture and utilize the important terms within an area. For example, the concept of "attainment level" can be organized as a multilevel taxonomy that is fully editable and extensible with the ability to select synonyms as needed. Figure 4 shows part of an attainment taxonomy within a design environment. It shows three key taxonomies often used in building learning outcomes.

**Figure 4.** Part of an "attainment" taxonomy for building ILOs.

Similar approaches can lead to taxonomies for all the elements within Figure 3. Of particular interest are the discipline-specific knowledge taxonomies. Figure 5 shows a taxonomy for Chemical Engineering that is both extensible and editable within the design environment. It consists of levels that express a *domain*, *sub-domain* and *topic* taxonomy. This is clearly a multilevel view.

**Figure 5.** Knowledge taxonomies for Chemical Engineering (domain/sub-domain/topic).

#### 3.2.2. Learning Unit Design

The key object in any curriculum design is the learning unit. This requires the integration of many elements in Figure 3:


This object is complex and could be regarded as analogous to a process engineering unit of a flowsheet. Learning outcomes can be built by using a defined but flexible syntax that uses the various taxonomies displayed in Figure 3. Structuring the ILOs provides significant power to reason over the completed curriculum and also visualize the characteristics of the curriculum. An example ILO from an advanced modeling course introducing the theory and practice of hybrid modeling is shown in Figure 6.

**Figure 6.** Example of intended learning outcomes (ILOs) (left) and the structured text ILO builder (right).

Moreover, the various tabs on the course profile show the other course-building options that are used in fully specifying the course characteristics. They include assessment items, learning activities, prior learning for the course and important mappings to professional competences.

#### *3.3. Pathways*

Curriculum designs ultimately lay down a series of learning pathways, professional attributes and competency development. As noted previously in the introduction, the program schema elements are developed incrementally across the length of the program. The representation and tracking of this incremental development are a crucial part of the design process.

For example, we often desire to understand the linkages within years and across years. This is similar to process streams passing through operating units. A simple visual illustration of course integration via learning outcomes is shown in Figure 7 [17].

**Figure 7.** Example of knowledge domain outcome linkages for process system courses. (By permission of Engineers Australia [17]).

#### **4. Curriculum Design:** *The Journey Maker* **Computer-Aided Design Environment**

In the previous sections, we described the importance and challenges of complex curriculum design. The key goals were discussed, as well as the design information and processes required to build curricula. The development of a comprehensive web-based design tool which is adaptable, extensible and editable provides for design across all disciplines across higher education. Importantly, it was also developed for student use in building personalized curricula and visualizing the impact of choices. This is particularly important in the Arts and Humanities disciplines.

The functionality of this curriculum environment includes the following [17]:


The majority of this functionality was developed in 2011/2012 as a standalone application [17]. The further development expanded activity across all university faculties, resulting in a web-based environment called *The Journey Maker*. The environment has two major components: a design component and a visualization component called "The Visual Journey". Figure 8 shows the entry page to the development environment, showing the design functions available for use in building and visualizing curricula.

*Processes* **2020.**, *8*, 373

**Figure 8.** Overview of *The Journey Maker* design-environment functionality.

Elements of this environment are illustrated in the next section, using the design of some selected learning units, as well as showing an overall five-year curriculum design.

#### **5. Case Study in Curriculum Design: A Five-Year Chemical Engineering Program**

As a case study to display the development of a curriculum, we use the design of a five-year combined Bachelor and Masters Chemical Engineering program (BE/ME) at The University of Queensland. This is sufficiently complex, to demonstrate the design tools and the visualization tools that give deep insight into the design. The overall program structure is shown in Figure 9. This program typically has four learning units per semester, including elective options, which are shown as "+" options in each semester. These are chosen by students to fulfill their degree requirements, build specializations and address personal interests. There are 10 semesters in this program, with industry placements in Year 4, Semester 2.

**Figure 9.** Compulsory course requirements for five-year BE/ME program.

As shown in Section 3.2.2, all the course profiles have been built by using the course builder, and so they have a rich description of the learning within each module and also the connectivity throughout the curriculum. It is possible to investigate this particular design through the use of the visualization facilities within "The Visual Journey" web environment.

#### *Overall Curriculum Characteristics*

The visualization is able to view the overall picture of the curriculum as seen through the ILOs. Figure 10 shows the relative frequency of ILOs across the principal domains of the program. Clearly, the knowledge domain of Chemical Engineering dominates, but we also see the relative emphases around professional competences and the basic sciences: mathematics and chemistry. These insights might lead to unit-level or curriculum-level redesigns, depending on overall outcomes requirements. The visualizations can be presented in many ways: frequency, weighted frequencies and the like.

**Figure 10.** ILO frequencies for major knowledge domains across whole curriculum.

For further insights, the individual domains can be expanded, as seen in Figure 11, which focuses on the sub-domains within the Chemical Engineering domain. The thumbnail of the whole curriculum can be viewed to see where a specific focus might occur. In this case, process design can be seen within the ILOs across a series of courses starting in Year 1.

**Figure 11.** Lower-level domain knowledge visualization and location within the curriculum.

This form of visualization, and many others, can also be applied to other aspects of the learning outcomes, such as cognitive objectives as described by Bloom's taxonomy. Figure 12 shows the distribution of ILOs where "synthesis" is a major objective within courses. In particular, the issue of synthesis described by "design" can be seen in later years of the curriculum. If some redistribution of that attainment to earlier years was considered important, then the visualization helps in redesign.

**Figure 12.** Attainment levels and distribution across the curriculum.

As well as considering issues around where knowledge, use of knowledge and personal attributes are developed, it is important to understand the linkages through the curriculum. This is particularly important as curriculum evolves, and in many cases, fragmentation occurs over time due to loss of integration. Figure 13 shows the connections of a particular ILO into following courses. Loss of that ILO due to changes in the learning unit could be important. The alluvial plots also allow unit instructors to see the position of the unit in relation to other units from a learning perspective. They see the required prior learning and also the future use of learning outcomes.

**Figure 13.** Tracking the linkages of a specific learning outcome through courses in the curriculum.

Figure 14 shows how course integration occurs across the curriculum. It shows key linkages along important learning pathways and identifies what might be "orphaned" courses within the curriculum.

**Figure 14.** Tracking linkages between courses that have specified required prior learning.

Many other tabular and graphical views can be generated due to the structured descriptions and taxonomies embedded within the design environment. Within the approach, there is a very extensive subsystem for designing the assessment strategies. Figure 3 shows the extent of all the curriculum components, which cannot all be discussed in this paper.

These provide the ability for any discipline to develop curricula that is sharable, easily updated and extensible for the specific issues often dealt with in higher education programs.

The effects of using these formalized approaches, as described in these developments, are that now staff and students can clearly see the interconnections within curricula and are able to see the flow of knowledge and skills that are needed for future learning. It also allows curriculum changes to be easily mapped and tracked over time, providing easy access to professional accreditation organizations to see the impact of those changes. The ability to visualize learning outcomes and the forms of assessment used to provide evidence of attainment levels is an important aspect of the design.

Due the structured, extensible nature of all the curriculum building blocks outlined in Figure 3, it is relatively easy to export a wide range of reports, visualizations, curriculum data and statistics for numerous purposes.

#### **6. Conclusions**

This work has shown the use of certain PSE principles in developing complex curricula. The design of curricula resembles in many ways the basic ideas of process design and flow sheet development, since curricula design is an educational process. Through the understanding of the crucial building blocks for curricula and the structured manner of developing learning units, it is possible to produce whole curricula designs that capture the many characteristics that make up complex learning environments.

By using structured information management approaches, the final designs can be visualized to understand the whole integrated curriculum. The fundamental concepts, organization and deployment started within Chemical Engineering, but they have now been adopted and used in many other disciplines, including nursing, science, agriculture, pharmacy, veterinary science, medicine and philosophy. It is significant that the idea of the importance of an integrative approach rather than disparate designs was at the heart of the PSE efforts of Roger Sargent from the very beginning of his long and distinguished career.

The application of PSE ideas to bring an integrative approach to curriculum design now yields similar benefits to many other disciplines that have already been realized within the PSE community.

**Author Contributions:** Conceptualization, I.C. and G.B.; methodology, I.C. and G.B.; software, eLIPSE.; formal analysis, I.C. and G.B.; investigation, I.C. and G.B.; resources, I.C. and G.B.; data curation, I.C.; writing—Original draft preparation, I.C., G.B.; writing—Review and editing, I.C. and G.B.; project administration, I.C. and G.B.; funding acquisition, I.C. and G.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** We acknowledge funding from The University of Queensland, Technology Enabled Learning (TEL) Grant scheme.

**Acknowledgments:** We acknowledge the web development work of the Centre for eLearning Innovations and Partnerships in Science and Engineering (eLIPSE). We thank Erzsébet Németh (E.N.) for her work on the visualization aspects of curricula.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Processes* Editorial Office E-mail: processes@mdpi.com www.mdpi.com/journal/processes

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18