Biobjective Optimization Model Considering Risk and Proﬁt for the Multienterprise Layout Design in Village-Level Industrial Parks in China

: With the advent and development of Industry 4.0 and 5.0, manufacturing modes have changed and numerous newly complicated and integrated village-level industrial parks have emerged in the Southeast of China, where several enterprises are gathered in the same multistory building. The number of ﬂoors and surrounding enterprises can have an impact on accident risk. To reduce the overall risk level of industrial parks, the layout of enterprises with different risks needs to be well designed and optimized. However, to date, limited studies have been conducted to emphatically consider safety and optimize the enterprise layout at an industrial area level, and most studies focus on the cost of the layout. Therefore, this study proposed three biobjective mathematical optimization models to obtain the trade-off between minimizing risk and maximizing rental proﬁt. Risk factors include the enterprise location and the association risk; the enterprise inherent safety risks are not considered. To solve this problem, a speciﬁc linearization strategy was proposed and an epsilon-constraint method was applied to obtain Pareto-optimal solutions. Subsequently, an industrial park in Shunde, China, was considered as a case study to verify the performance of the proposed models and methods. Finally, a sensitivity analysis of critical parameters was conducted. The critical factors inﬂuencing the objective functions were also analyzed to provide valuable managerial insights.


Introduction
Under the background of China's reform and opening up, many "workshop" villagelevel industrial parks have appeared in the Pearl River Delta region of China.Historically, village-level industrial parks have contributed significantly to the local economy.However, with the development of society, especially with the advent of Industry 4.0 [1] and "Made in China 2025" [2,3], these types of village-level industrial parks, with low buildings, unreasonable layout, and a low land utilization rate, has not met the needs of social development.Consequently, a new type of intensive village industrial park arises, which consists of multiple multistory standard buildings.In the new village-level industrial parks, different types of enterprises are concentrated, and are even located in the same building (see Figure 1).When an accident occurs in an enterprise, it may cause harm to neighboring enterprises and employees [4], or even cause a domino effect.Accordingly, how to reduce the risk level of new-type village industrial parks by the scientific layout of enterprises has become the focus of research.Our study optimizes the layout from a park management perspective.In addition, the cost invested by the park management is concentrated in the construction stage of the park, and the profit mainly comes from rent affected by the layout results.Therefore, it is necessary to consider rental profit.
models.The optimization of three-dimensional layouts at the level of industrial parks, with enterprises as the smallest layout unit, is, however, a topic that has received little research.In addition, most studies have focused on layout costs [9], such as pipeline and land costs.Nevertheless, rental profit has rarely been considered in previous studies.Furthermore, with the development of society, safety problems have become a critical issue in supporting the sustainable development of society.Industry 5.0, proposed by the European Union, also emphasizes people-oriented and sustainable development [10].Although there have been significant advances in the field, there are still research gaps in current research.First, there is a lack of research on the layout problem considering multibuilding, multistory, and multienterprise factors.Secondly, there is little research on the trade-off between risk and rental profits.To address the above gaps, this study used three biobjective mathematical optimization models to obtain the trade-off between minimizing risk and maximizing rental profit.In this sense, this study mainly focuses on the following key questions: (1) What type of risk needs to be considered?Do different types of risks have a significant impact on layout results?(2) How can we formulate a biobjective model for the three-dimensional layout problem considering multibuilding, multistory, and multienterprise factors?(3) How can we linearize and solve the proposed nonlinear mathematical programming model?(4) What are the managerial insights for decision-makers in handling this new type of layout problem?
To solve the problem of enterprise layout optimization with the aforementioned questions, mathematical programming modeling was performed in this study.Compared with existing studies, this study makes the following main contributions: (a) it focuses on the two objectives of minimizing risk and maximizing rental profits; (b) it considers a The layout problem has a long history, and several studies have been conducted on the layout of facilities inside a plant [5] presented a combined quantitative and qualitative approach to the facilities layout problem.A mixed integer nonlinear programming (MINLP) model was created by Penteado and Ciric [6] to determine the ideal layout of process equipment in a chemical plant.However, the issue of the multistory layout was not taken into account in these investigations.Research on multistory layouts has gained popularity as the number of multistory workshops has increased.In the context of the facility layout problem, Kia et al. [7], and Latifi et al. [8], offered two three-dimensional models.The optimization of three-dimensional layouts at the level of industrial parks, with enterprises as the smallest layout unit, is, however, a topic that has received little research.In addition, most studies have focused on layout costs [9], such as pipeline and land costs.Nevertheless, rental profit has rarely been considered in previous studies.Furthermore, with the development of society, safety problems have become a critical issue in supporting the sustainable development of society.Industry 5.0, proposed by the European Union, also emphasizes people-oriented and sustainable development [10].
Although there have been significant advances in the field, there are still research gaps in current research.First, there is a lack of research on the layout problem considering multibuilding, multistory, and multienterprise factors.Secondly, there is little research on the trade-off between risk and rental profits.To address the above gaps, this study used three biobjective mathematical optimization models to obtain the trade-off between minimizing risk and maximizing rental profit.In this sense, this study mainly focuses on the following key questions: (1) What type of risk needs to be considered?Do different types of risks have a significant impact on layout results?(2) How can we formulate a biobjective model for the three-dimensional layout problem considering multibuilding, multistory, and multienterprise factors?(3) How can we linearize and solve the proposed nonlinear mathematical programming model?(4) What are the managerial insights for decision-makers in handling this new type of layout problem?
To solve the problem of enterprise layout optimization with the aforementioned questions, mathematical programming modeling was performed in this study.Compared with existing studies, this study makes the following main contributions: (a) it focuses on the two objectives of minimizing risk and maximizing rental profits; (b) it considers a three-dimensional layout optimization problem with multibuilding, multistory, and multienterprise factors, where the calculation scale of this study is larger than that of the two-dimensional layout design cases; and (c) it considers the risks of the park from multiple perspectives, such as the risk of enterprises on different floors and the mutual influence among enterprises.The remainder of this paper is structured as follows.The key issues of this study are reviewed in Section 2, where research gaps are also discussed.The enterprise layout problem, with its objectives and presumptions, is described in Section 3. To address this problem, biobjective mathematical models are constructed in Section 4, and the linearized and solution strategies are developed in Section 5.In a case study presented in Section 6, the computed findings support the model's efficacy.Section 7 summarizes the results and highlights future directions.

Layout Object
Layout design optimization impacts the design of many products and industrial systems [11].The research objects were mainly facilities, products, and enterprises/factories [12].Existing studies mostly focus on the facility layout problem [12,13].Patsiatzis et al. [14] studied the single-floor process plant layout problem, considering the piping cost, property damage, and cost of protection devices.Xu et al. [15] and Medina-Herrera et al. [16] proposed mathematical models for facility layout in chemical plants.These studies take facilities as the smallest unit of the layout.Our study aims to optimize the layout of different enterprises in an industrial park, that is, the overall regional layout.Wang et al. [17] studied the effects of multiple hazard sources on plant layout in industrial parks.Wu and Wang [18] proposed a systematic optimization method for the chemical area layout, considering safety and environmental issues.Traditional regional layouts are generally single-floor and do not involve multistory buildings.However, enterprises are becoming increasingly intensive.In such a situation, as a secondary allocation problem, the enterprise layout design problem represents an important combinatorial optimization problem, which is an urgent task to be solved.

The Multistory Layout Problem
Many studies have been conducted on two-dimensional facility layout problems [19,20].With an increase in new industrial parks and high-rise buildings, research on multistory layout problems has gradually become a trend.A multistory layout needs to consider the impact of different floors on the layout objective.Fewer studies have investigated the layout of multistory plant facilities [7].For instance, O'Neill et al. [21] introduced the z-axis on a two-dimensional basis for the facilities layout inside a factory, which inspired this work.Latifi et al. [8] used positive variables to describe the position of the center of the facility.The layout location of the park is fixed.This study chooses binary variables to represent the locations of the enterprises and constructs multiple mathematical models in the form of a three-dimensional layout to study the overall layout of enterprises in industrial parks.

The Layout Principle
Existing industrial layout problems focus on the economy, and aim at minimizing the total layout cost [20,22].Following on, Rw et al. [23] set the pipeline cost as an objective function to optimize the layout.Earlier et al. [24] considered reducing material handling costs among all departments, and the cost of rent.In 2021, Hosseini et al. [25].used material handling costs and machine rearrangement costs to determine the best layout.These studies conduct layout optimization from the perspective of enterprises.Similarly, minimizing pipeline investment costs, pump power costs, land costs, and surface construction costs, were considered as the goals of optimal layout [26].Our study takes the enterprise as the smallest layout unit and considers the interests of the park management to optimize the layout.The cost invested by the park management is concentrated in the construction stage of the park, and the profit mainly comes from rent.Unlike previous studies that considered minimizing costs, this study sets maximizing rental profits as an objective function.A similar study exists that focuses on profit [27].
Previous studies involving safety have considered risk factors such as cost or layout constraints [28].Further, Vazquez-Roman et al. [29] used land, pipeline, and risk costs as objective functions, and also considered the effects of toxic release constraints.With the development of society, the importance of safety is increasing as it is a critical issue to support sustainable development.Safety and risk assessment, which are presently in-creasing, were seldom considered in previous studies in the past century [30].Interestingly, Wang et al. [31] proposed a multiobjective optimization method to achieve different tradeoffs between the economy and security.This study concerns safety issues and rental profit in the layout process.

The Mathematical Model
Mathematical programming methods have been widely used to solve industrial layout problems [32].For instance, de Lira-Flores et al. [33] proposed a MINLP method to obtain the optimal facility layout for processing plants.Further, Brunoro Ahumada et al. [34] incorporated the risk diagram and safety distance into a mixed-integer linear programming (MILP) model for facility layout optimization.None of the studies considered time.Many sustainable dynamic layout models have also been developed in recent years [35][36][37].Generally, the location of an enterprise in an industry park does not change frequently; therefore, our models did not add time variables.For the linear programming model, the branch and bound and simplex algorithms can be used for accurate solutions [38].The nonlinear programming model can be transformed into a linear programming model for an accurate solution, and the metaheuristic algorithm can also be used to provide solutions [39,40].For instance, Wang et al. [41] optimized the layout of industrial parks based on a genetic algorithm.Further, Alves et al. [42] used simulated annealing to move units throughout the chemical plant until a workable layout was found to minimize the impact on the surrounding community, in the event of an accident.The optimal solution of the mathematical programming model obtained using the metaheuristic algorithm is the approximate optimal solution.Our study proposes a specific linearization strategy to use the exact solution method.

Summary
In addition to the above discussion, Table 1 summarizes a review of relevant studies.Because this study focuses on the safety and rental profit in park layout, international journals that have the highest relevant contributions have been researched.These include: Safety Science; Process Safety and Environmental Protection; the Annals of Operations Research; the Journal of Loss Prevention in the Process Industries; and Computers and Chemical Engineering.As can be seen from Table 1, the vast majority of studies consider costs; none focus on rental profits.In addition, numerous studies are devoted to optimizing the two-dimensional layout of the facility.Therefore, our study focuses on a three-dimensional layout optimization problem, with multibuilding, multistory, and multienterprise factors.Three biobjective models are presented, considering the risks from multiple perspectives and rental profit.To solve this problem, a linearization strategy is proposed, and the epsilon-constraint method is adopted to solve the reformulated models.This new type of village industrial park has several intensive modern standard plants.There are multiple multistory buildings in the industrial park.They are separately leased to multiple independent enterprises, and unified management is implemented by the park property.Enterprises are densely distributed in industrial parks, which significantly improves land utilization rates.There are different types of enterprises in the park.For instance, common industries in Shunde Industrial Park include: machinery and equipment; household appliances; packaging and printing; electronic information; biological medicine; and other industries.Each enterprise has different sources of hazard and different safety risk values, located on different floors.Meanwhile, enterprises located in the same building share the building infrastructure, are close to each other, and influence each other.In general, the determining factors of the overall safety risk of the industrial park include: (a) the enterprise inherent safety risk; (b) the enterprise location safety risk; (c) the enterprise association safety risk; and (d) the domino effect.

Enterprise Inherent Safety Risks
Generally, different types of enterprise are located in the industrial park, where many types of accidents may occur, such as fire, explosion, object strike, mechanical injury, poisoning, and asphyxiation.The risk points within the enterprise should be investigated and identified and the type of accident, based on relevant standards, should be determined.An appropriate risk assessment method was adopted to determine the security risk value of an enterprise.The value of the enterprise's inherent safety risk only considers the risk factors of the enterprise; the influence of surrounding enterprises is not considered and does not vary with the location change.In addition, owing to the different nature of the industry, enterprises generally have different risk factors, and their inherent safety risk values differ.

Enterprise Location Safety Risks
From the qualitative perspective of engineering experience, the risk values of enterprises located on different floors are different.The quantitative evaluation method can be used to evaluate the enterprise location safety risk.The objective of this method is to evaluate the possibility of potential accidents and the severity of the consequences; thus, the risk value is the combination of the above likelihood and severity.A change in the location of the same enterprise may result in a change in the likelihood or severity of an accident.The enterprise location safety risk can be obtained by the following formula: where K As is the risk value of Enterprise A on floor s, L As is the possibility of Enterprise A having an accident on the floor s, and C As is the loss caused by the accident of Enterprise A on floor s.
This study provides some visual examples to better illustrate the location safety risk.Enterprises with heavy equipment are almost impossible to collapse when they are located on the ground floor of a building; however, there is a risk of collapse when they are located at the top of the building.Enterprises with an explosion risk should be arranged on the top floor of the building, as the spread range and consequences caused by an explosion accident are more likely to be smaller than those on the bottom floor, as shown in Figure 2. In this study, the difference in the risk values of enterprises located on different floors is considered, that is, the location risk of enterprises is considered.
where  is the risk value of Enterprise A on floor ,  is the possi A having an accident on the floor , and  is the loss caused by the prise A on floor .
This study provides some visual examples to better illustrate the l Enterprises with heavy equipment are almost impossible to collapse cated on the ground floor of a building; however, there is a risk of colla located at the top of the building.Enterprises with an explosion risk s on the top floor of the building, as the spread range and consequence plosion accident are more likely to be smaller than those on the bottom Figure 2. In this study, the difference in the risk values of enterprises l floors is considered, that is, the location risk of enterprises is considere

Enterprise Association Safety Risks
Enterprises located in the same building in the park may influence each other because they are close to each other and share buildings.Enterprises that produce or use corrosive solutions and gases will increase the risk of accidents if they are surrounded by enterprises that emit large amounts of dust.Enterprises with a potential gas leak will bring additional risks to surrounding enterprises.An explosion in one enterprise may cause an explosion within another potentially explosive enterprise.The enterprise association safety risks can be obtained by the following formula [46]: where R AB is the association risk of Enterprise A to Enterprise B, P A is the possibility that Enterprise A will have an accident, P AB represents the likelihood that an event in Enterprise A will result in a second accident in Enterprise B, and C B indicates the damages caused by the second accident in Enterprise B.
In general, as shown in Figure 3, the influence between enterprises can be roughly divided into three categories: (a) the symmetrical interaction relationship: Enterprises A and B have the same association safety risks; (b) the unilateral relationship: the impact of Enterprise B on Enterprise A is negligible; and (c) the asymmetric interaction relationship: Enterprises A and B have different association safety risks.There are often many enterprises in industrial parks, and the association safety risk among enterprises significantly influences the overall risk level of the park.The width of the arrows represents different degrees of influence.
divided into three categories: (a) the symmetrical interaction relationship: Enterprises A and B have the same association safety risks; (b) the unilateral relationship: the impact of Enterprise B on Enterprise A is negligible; and (c) the asymmetric interaction relationship: Enterprises A and B have different association safety risks.There are often many enterprises in industrial parks, and the association safety risk among enterprises significantly influences the overall risk level of the park.The width of the arrows represents different degrees of influence.

Domino Effect
The domino effect [47,48] is the idea that in an interconnected system a small initial energy can lead to a series of chain reactions.The geographical closeness between enterprises in industrial parks, especially among enterprises with large-hazard sources, leads to the enhancement of the domino effect.That is, when an accident occurs in one enterprise, it often directly or indirectly influences one or more enterprises and compromises the safety of other enterprises.For example, an accident in an enterprise leads to the accidental release of energy, such as fire, explosion, or gas leakage, which leads to the failure of facilities in other nearby enterprises, resulting in a series of accidents.The consequences of the accident are larger, and the loss is more serious, than that of the initial accident.This chain reaction increases the severity of accidents.The domino effect is a key issue to be

Domino Effect
The domino effect [47,48] is the idea that in an interconnected system a small initial energy can lead to a series of chain reactions.The geographical closeness between enterprises in industrial parks, especially among enterprises with large-hazard sources, leads to the enhancement of the domino effect.That is, when an accident occurs in one enterprise, it often directly or indirectly influences one or more enterprises and compromises the safety of other enterprises.For example, an accident in an enterprise leads to the accidental release of energy, such as fire, explosion, or gas leakage, which leads to the failure of facilities in other nearby enterprises, resulting in a series of accidents.The consequences of the accident are larger, and the loss is more serious, than that of the initial accident.This chain reaction increases the severity of accidents.The domino effect is a key issue to be considered when enterprises optimize the layout of industrial parks.The strength of the domino effect can be represented by the enterprise association safety risk.
The enterprises' inherent risks generally do not change significantly with the change of location, but the layout of enterprises in the park will have a significant impact on the location risks and association risks.Therefore, our study considers the location and association risks during the layout process.

Profit
In addition to safety issues, profit is also the focus of park management, which is affected by many factors, such as factory rent, management costs, and government support.Among these, rent is affected by the layout results and is an important source of profit for the park.The new-type village industrial park generally consists of several multistory buildings; different enterprises have different affordability values for each floor.For example, enterprises with a larger flow of goods are more willing to pay higher rent for the first floor; therefore, it is possible to maximize rental profit by optimizing the park layout.

Objectives
In general, the optimized layout of enterprises aims to obtain the trade-off between minimizing the overall risk and maximizing the rental profit.Based on different types of risk, we have established three biobjective mathematical programming models: (a) Model 1: the first goal aims to achieve the minimum enterprise location safety risk and the second goal is to maximize the rental profit; (b) Model 2: the first goal aims to achieve the minimum enterprise association safety risk and the second goal is to maximize the rental profit; (c) Model 3: the first goal aims to achieve the minimum enterprise location and association safety risk and the second goal is to maximize the rental profit.

Presumptions
Before presenting the proposed biobjective mathematical models, four assumptions are postulated: (1) Enterprises are treated as plane figures of different sizes, which can be placed in the free spaces of any building in the industrial park.(2) An enterprise can only be located on one floor of the building because small-and medium-sized enterprises almost exclusively occupy village-level industrial parks.
The few enterprises that need to be located on multiple floors are treated as multiple independent enterprises.(3) The situation in which the location of some enterprises are determined in advance is allowed.Additionally, the impact of these enterprises on the overall risk value and rental profit of the park is considered.Enterprises need to be entirely housed in an industrial park and do not overlap with each other.
The association risk can only exist between two enterprises if they are housed in the same building.The influences between enterprises situated in different buildings were not considered in this study; we can find a similar presumption in Liu et al. [27].

The Safety-Oriented Optimization Model
This section provides three biobjective mathematical optimization models to formulate three problems with different concerns, in the following three subsections.The notations used in these three models are first provided before presenting the three models.Please note that some notations are the same for the three models.

. Model Formulation
The model 1 is defined using Equations ( 3)-( 6). Min The first Objective Function (3) is to minimize the total location safety risk value, and the second Objective Function ( 4) is to maximize the total rental profit.Constraint (5) guarantees that the sum of the enterprise area of each floor cannot exceed this monolayer area.Constraint (6) requires an enterprise to have only one location and cannot be located in more than one location.It is worth noting that when the location of some enterprises has been determined, it is necessary to change Constraint (5); that is, the layer area to the space area that can be laid out.A ij → Associated risk value with enterprises i and j.Decision variables: 1 Enterprises i and j in building b 0 otherwise

Model Formulation
The model 2 is defined using Equations ( 3)-( 6). Min Max The first Objective Function ( 7) is to minimize the overall associated safety risk value, and the second Objective Function ( 8) is the maximization of the total rental profit.Constraint (9) guarantees that the sum of the enterprise area of each floor cannot exceed this monolayer area.Constraint (10) and (11) restrict each enterprise from being placed in one location.Constraint (12) guarantees that each enterprise can only be assigned to one building.To explain Constraint (13), a clear illustration is provided.All the enterprises are assumed to be in the same building.Thus, the enterprise association risk matrix (i.e., Equation ( 14)) can be obtained, where A ij represents the associated safety risk caused by the influence of enterprise i on enterprise j.The symmetrical elements (i.e., A ij and A ji ) do not have to be equal.
The possibility of association safety risk (i.e., A ij and A ji ) exists only if enterprises i and j are located in the same building.Accordingly, x bij was selected as the decision variable determined by y bsi and z bsj .That is, x bij can be 1 only when both enterprises i and j exist in building b.The relationships among decision variables (i.e., y bsi , z bsj , and x bij ) can be divided into four types, as summarized in Table 2.
The first Objective Function (15) is to minimize the sum of the location safety and association safety risk value.

Linearization and Solution Methodologies
In this section, the proposed biobjective models can be effectively solved using the IBM CPLEX Optimizer, due to the development of specific linearization and solution approaches.

Linearization Strategy
Because Constraint (13) in the proposed Models 2 and 3 is nonlinear, it is necessary to develop a linearization method.Although the heuristic method is significant in solving nonlinear problems, the difficulty in searching for the global optimum is a significant challenge.Consequently, this study develops a specific linearization method that introduces several critical parameters into the model to reformulate the proposed Models 2 and 3. To linearize the proposed nonlinear Models 2 and 3, the key problem is the removal of the MIN function in Constraint (13).In this study, a large positive number M and a natural number ρ are introduced into the proposed Models 2 and 3 such that the nonlinear Constraint (13) can be substituted by two linear constraints, which are given by: where ρ∈(1, 2) and M is a large number.Therefore, the proposed Models 2 and 3 can be reformulated as linear mathematical optimization models.Before presenting the re-formulated linear mathematical optimization models, the validity of the substitute linear constraints is provided as follows: Theorem 1. Constraint ( 13) is equivalent to the Inequalities in ( 16) and ( 17), while ρ∈ (1,2).
Proof.As presented in Table 2 and Constraint (13), x bij is the minimizer value between ∑ s∈S y bsi and ∑ s∈S z bsj .Instead of using the MIN function in Constraint ( 13), a combination of Inequalities ( 16) and ( 17) is used.It suffices to show x bij assumes the minimizer value of ∑ s∈S y bsi and ∑ s∈S z bsj whereas ρ∈(1, 2).Because ∑ s∈S y bsi and ∑ s∈S z bsj can only assume binary values (i.e., 0 or 1), as shown in Table 2, four different combinations exist (i.e., Types 1-4).Thus, we validated the effectiveness of the combination of Inequalities in ( 16) and ( 17) in these four combinations.
(1) Type 1 Given ∑ s∈S y bsi = 1 and ∑ s∈S z bsj = 1, we have the below inequalities: Because ρ ∈ (1, 2), it is easy to see that 2−ρ M is a small positive value.Let us use ε + notation to represent 2−ρ M and then we have ε + ∈ (0, 1).Please note that M is a large positive value.Since the above two inequalities need to be satisfied, then we have: Since x bij can only assume binary values, x bij is 1 with the above restriction in (20).
(2) Types 2 and 3 Given ∑ s∈S y bsi = 1 and ∑ s∈S z bsj = 0 OR ∑ s∈S y bsi = 0 and ∑ s∈S z bsj = 1, we have the below inequalities: Because ρ∈(1, 2), it is easy to see that 1−ρ M is a negative value.Let us use notation ε − to represent 1−ρ  M and then we have ε − ∈ (−1, 0).Since the above two inequalities need to be simultaneously satisfied, then we have: Since x bij can only assume binary values, x bij is 0 with the above restriction in (23).
(3) Type 4 Given ∑ s∈S y bsi = 0 and ∑ s∈S z bsj = 0, we have the below inequalities: Because ρ∈(1, 2), it is easy to see that − ρ M is a small positive value.Let us use ε * notation to represent − ρ M and then we have ε * ∈ (−1, 0).Since the above two inequalities need to be satisfied, then we have: Since x bij can only assume binary values, x bij is 0 with the above restriction in (26).Then, we complete the proof.Theorem 2. To ensure that the combination in Inequalities ( 16) and ( 17) is equivalent to Constraint (13), the minimum value for Mis 2.
Proof.As stated in the above-mentioned proof for Theorem 1, it is assumed that M has a large positive value.To satisfy Inequalities ( 20), (23), and (26), it suffices to show that: Then, we have a combination of inequalities, which is given as below: Since ρ∈(1, 2), it is easy to see that M ≥ 2.Then, we complete the proof that the minimum value for M is 2.

The Epsilon-Constraint Method
The models proposed in this paper are three biobjective models.Epsilon-constraint, weighted sum, weighted metric, and lexicographic goal programming approaches are some of the traditional biobjective optimization methods that are commonly applied in practice.The weighted sum and weighted metric methods require the introduction of new parameters.The lexicographic goal programming approach is used when there exists a clear priority ordering amongst the goals to be achieved.The two goals we focus on have no obvious priority.The epsilon-constraint method is used in this work to provide the nondominated solutions because of the advantages, which can be found in Gao and Cao [49].The main idea of this method is to convert an objective function into a constraint and establish a connection between the two objectives by gradually decreasing or increasing the E 1 value [50].The three modes can be formulated as follows: Constraints ( 5) and ( 6) Constraints ( 9)-( 12), ( 16) and ( 17) Constraints ( 9)-( 12), ( 16) and ( 17) Taking Model 1 as an example, the steps of the epsilon-constraint method are as follows: (i) the minimum value of the single objective function Ψ 1 is calculated under Constraints ( 5) and ( 6); (ii) solve the Model 1 given in Constraints ( 5) and ( 6); (iii) substituting the layout result into Function (3) to obtain the corresponding nondominated solutions; and (iv) change the value of E 1 , and then repeat steps (ii) and (iii).Thus, the Pareto frontier is obtained.

Case Study
This study focuses on the layout of multiple enterprises in a village-level industrial park, and it is unique to consider multistory buildings.The effectiveness of the proposed models and methods was tested using a case study of a Shunde industrial park, and the following results are presented.This section comprises four subsections.In the first part, case data are introduced and processed.Subsequently, data are substituted into the three proposed models for calculation, and the results are presented.In part three, the sensitivity analysis was performed by changing the parameter values, and some management opinions were proposed based on the simulation results.The IBM ILOG CPLEX Optimization Studio (Version:12.6)was applied to solve the proposed biobjective models.All the experiments were run on a computer with an 11th Gen Intel(R) Core (TM) i5-1135G7 @ 2.40 GHz and 16 GB memory, with a Windows 10 operating system (Lenovo (Beijing) Co., Ltd., Beijing, China).

Data Introduction
In this subsection, the layout of an industrial park in Shunde is considered as a case study to evaluate the three proposed models.The building information for this industrial park is presented in Table 3.The industrial park has four buildings with five floors each, and 20 enterprises need to be rationally arranged in the park.The type of enterprise include: packaging and printing; machinery and equipment; household appliances; electronic information; biological medicine; and e-commerce.The area and number of floors of all enterprises are known.Given that it is challenging to obtain comprehensive and accurate data, we assumed the risk and rent parameters based on the type of settled enterprise, the local environment, and the local economy.

Illustrative Results
Here, the main results for the three-dimensional layout problem of the enterprise are presented by considering different risk types and rental profit.The value of E 1 increases by 50 each time, and the Pareto frontiers are shown in Figure 4.These two objective functions are generally in conflict with one another, as is evident from Figure 4. We cannot guarantee that the lowest level of risk and the highest rental profit will be achieved at the same time.The second objective function value worsens as the first objective function value is decreased.The new industrial park layout needs to pay attention to risks, but the park management will consider rental profit in the actual layout process.Therefore, the Pareto solution can provide a reference for decision-makers, who can choose the corresponding layout scheme according to the actual needs of the park.None of the nondominated solutions in the set can be stated to be preferable to others.Therefore, for the decisionmakers to determine the "most desired" alternative, more preference data is required.
The minimum overall risk levels of the park calculated using Models 1, 2, and 3 were 194, 1031, and 1225, respectively.The corresponding rental profits were CNY 1,544,125, 1,575,281, and 1,541,801.The layout results of the three scenarios are shown in Figure 5a,c,e, where I1-I20 represent 20 enterprises.Model 2 considers only the association safety risk, which only exists in enterprises in the same building.The association safety risk was unrelated to the floor number, so the lowest risk scenario is more likely to obtain a larger rental profit than the other models.The maximum rental profit calculated through the three models is the same (i.e., CNY 1,588,165).The specific layout results with the greatest rental profit are shown in Figure 5b,d,f.It can be found that in the case of a certain rental profit, there are multiple layout results.The layout results of the three models are not the same.

Sensitivity Analysis
Sensitivity analysis was performed by changing the key parameters in this subsection.In this analysis, the case data were used as reference values.The specific enterprise layout result obtained from the three models.

Sensitivity to Floors
The impact of the number of floors on the overall risk value, and the rental profit of the park, were analyzed by changing the number of floors from five to eight.The optimal Pareto solutions of three models with different floors are partially shown in Tables 4-6.It is obvious that the minimum risk value and maximum rental profit do not improve in any way with the different number of floors.Model 1 only considers the location risk, which is affected by the number of floors.The minimum risk value has not changed because there is enough layout space.The location safety risk is not considered in Model 2, and the floor constraint in Model 2 only reflects the area constraint.With sufficient layout space, the risk value in Model 2 is almost independent of the number of floors.It is foreseeable that increasing the number of floors has no impact on maximum rental profit because enterprises are more willing to pay higher rents for the lower floors.The comparisons of the Pareto fronts of the three models are shown in Figure 6.In Model 2, the four Pareto frontiers almost overlap.Furthermore, it is interesting to find that increasing the number of floors with a fixed risk value has the potential to have a negative impact on rental profits.

Sensitivity Analysis
Sensitivity analysis was performed by changing the key parameters in this subsection.In this analysis, the case data were used as reference values.

Sensitivity to Floors
The impact of the number of floors on the overall risk value, and the rental profit of the park, were analyzed by changing the number of floors from five to eight.The optimal Pareto solutions of three models with different floors are partially shown in Tables 4-6.It is obvious that the minimum risk value and maximum rental profit do not improve in any way with the different number of floors.Model 1 only considers the location risk, which is affected by the number of floors.The minimum risk value has not changed because there is enough layout space.The location safety risk is not considered in Model 2, and the floor constraint in Model 2 only reflects the area constraint.With sufficient layout space, the risk value in Model 2 is almost independent of the number of floors.It is foreseeable that increasing the number of floors has no impact on maximum rental profit because enterprises are more willing to pay higher rents for the lower floors.
The comparisons of the Pareto fronts of the three models are shown in Figure 6.In Model 2, the four Pareto frontiers almost overlap.Furthermore, it is interesting to find that increasing the number of floors with a fixed risk value has the potential to have a negative impact on rental profits.

Sensitivity to the Number of Buildings
In different parks there were differences in the number of buildings.Based on the case data, one building was reduced and one building was added, to analyze the change in the overall risk value and profit of the park.The Pareto frontiers are shown in Figure 7, and the comparisons of the Pareto fronts of the three models are shown in Figure 8. Model 1 is the least sensitive to changes in the number of buildings.In Model 1, the minimum value of the location risk remains unchanged for different building numbers, indicating that the floor space is relatively sufficient when the number of buildings is three, four, and five.As the number of buildings increases, the minimum risk values calculated by Models 2 and 3 are significantly reduced.Only enterprises located in the same building may be associated with safety risks.Therefore, the number of buildings in the park significantly influences the association safety risk value.We can assume that there are sufficient buildings in extreme cases, such that each enterprise can be independently distributed in different buildings.In this case, the association safety risk value of enterprises in the park was 0. As the number of buildings increases, the maximum rental profit calculated by all three models rises by the same amount.influences the association safety risk value.We can assume that there are sufficient buildings in extreme cases, such that each enterprise can be independently distributed in different buildings.In this case, the association safety risk value of enterprises in the park was 0. As the number of buildings increases, the maximum rental profit calculated by all three models rises by the same amount.
Enterprises prefer to lease the lower floors, and the increase in the number of buildings provides more layout space for higher rental profit.

Sensitivity to the Number of Enterprises
The number of enterprises is changed from 18 to 26, and the resulting risk value and rental profit are output.The Pareto frontiers for the three models are presented in Figure 9, and the variations of the Pareto frontiers with a different number of enterprises are shown in Figure 10.It is evident from the figures that as the number of enterprises increase, the minimum risk value and maximum rental profit increases.For example, when the number of enterprises increases from 18 to 26, the minimum risk value calculated using Model 2 increases from 661 to 2297; the maximum rental profit calculated using Model 2 increases from CNY 1,399,245 to CNY 1,734,798.This situation is easy to explain.The increase in the number of enterprises will lead to a more intensive distribution of enterprises, and the choice of space of each enterprise will be reduced such that the total risk value may increase.Moreover, the enterprise itself has a location risk and associated risk with other enterprises.Adding enterprises involves adding new risks and rental profits to the park.
In the calculation process of Model 1, we observed an interesting phenomenon: the increased overall risk value is the lowest location safety risk of the additional enterprises.Therefore, for Model 1, the layout space is relatively affluent.Models 2 and 3 are sensitive to changes in the enterprise.This also indicates that an increase in the number of enterprises still dramatically impacts the total associated risk value in the case of sufficient layout space.Maximum rental profits do not show this phenomenon.Enterprises prefer to lease the lower floors, and the increase in the number of buildings provides more layout space for higher rental profit.

Sensitivity to the Number of Enterprises
The number of enterprises is changed from 18 to 26, and the resulting risk value and rental profit are output.The Pareto frontiers for the three models are presented in Figure 9, and the variations of the Pareto frontiers with a different number of enterprises are shown in Figure 10.It is evident from the figures that as the number of enterprises increase, the minimum risk value and maximum rental profit increases.For example, when the number of enterprises increases from 18 to 26, the minimum risk value calculated using Model 2 increases from 661 to 2297; the maximum rental profit calculated using Model 2 increases from CNY 1,399,245 to CNY 1,734,798.This situation is easy to explain.The increase in the number of enterprises will lead to a more intensive distribution of enterprises, and the choice of space of each enterprise will be reduced such that the total risk value may increase.Moreover, the enterprise itself has a location risk and associated risk with other enterprises.Adding enterprises involves adding new risks and rental profits to the park.
In the calculation process of Model 1, we observed an interesting phenomenon: the increased overall risk value is the lowest location safety risk of the additional enterprises.Therefore, for Model 1, the layout space is relatively affluent.Models 2 and 3 are sensitive to changes in the enterprise.This also indicates that an increase in the number of enterprises still dramatically impacts the total associated risk value in the case of sufficient layout space.Maximum rental profits do not show this phenomenon.

Sensitivity to Risk
The risk of enterprises entering the park is full of uncertainty.This study assumes an enterprise as an example to observe the change in enterprise layout, by changing its location risk and the association risk, when the overall risk of the park is minimal.The risk values for the parks mentioned in this subsection are minimized.
(1) Change in the enterprise location safety risk Enterprise 1 is used as an example to demonstrate the influence of a change in location safety risk on enterprise layout.The result of the enterprise layout in the park is calculated using Models 1 and 3, when the overall risk of the park is minimal.The location safety risk of Enterprise 1 is changed from 30, 30, 20, 1, 1 to 1, 1, 20, 30, and 30, respectively.The results of the layout with the minimum overall risk levels of the park are shown in Figures 11 and 12, and the minimum risk value of the park remains the original value.

Sensitivity to Risk
The risk of enterprises entering the park is full of uncertainty.This study assumes an enterprise as an example to observe the change in enterprise layout, by changing its location risk and the association risk, when the overall risk of the park is minimal.The risk values for the parks mentioned in this subsection are minimized.
(1) Change in the enterprise location safety risk Enterprise 1 is used as an example to demonstrate the influence of a change in location safety risk on enterprise layout.The result of the enterprise layout in the park is calculated using Models 1 and 3, when the overall risk of the park is minimal.The location safety risk of Enterprise 1 is changed from 30, 30, 20, 1, 1 to 1, 1, 20, 30, and 30, respectively.The results of the layout with the minimum overall risk levels of the park are shown in Figures 11 and 12, and the minimum risk value of the park remains the original value.
A comparison of Figures 5a and 11 show that with a change in the location risk of Enterprise 1, the original layout of Enterprise 1 changes from the 5th floor to the 1st floor.In Model 3, Enterprise 5 was changed from the original layout from the 5th floor to the 1st floor (see Figures 5e and 12).Moreover, the contribution of Enterprise 1 to the overall risk value of the park remains unchanged.It also shows that this model can optimize the layout of the enterprise according to the location safety risk to achieve the lowest overall risk value.(2) Change in the enterprise association safety risk According to Figure 5c,e, Enterprises 1 and 20 in Models 2 and 3 are arranged in the same building.To display the change in enterprise layout, the association safety risk value

Sensitivity to Risk
The risk of enterprises entering the park is full of uncertainty.This study assumes an enterprise as an example to observe the change in enterprise layout, by changing its location risk and the association risk, when the overall risk of the park is minimal.The risk values for the parks mentioned in this subsection are minimized.
(1) Change in the enterprise location safety risk Enterprise 1 is used as an example to demonstrate the influence of a change in location safety risk on enterprise layout.The result of the enterprise layout in the park is calculated using Models 1 and 3, when the overall risk of the park is minimal.The location safety risk of Enterprise 1 is changed from 30, 30, 20, 1, 1 to 1, 1, 20, 30, and 30, respectively.The results of the layout with the minimum overall risk levels of the park are shown in Figures 11 and 12, and the minimum risk value of the park remains the original value.
A comparison of Figures 5a and 11 show that with a change in the location risk of Enterprise 1, the original layout of Enterprise 1 changes from the 5th floor to the 1st floor.In Model 3, Enterprise 5 was changed from the original layout from the 5th floor to the 1st floor (see Figures 5e and 12).Moreover, the contribution of Enterprise 1 to the overall risk value of the park remains unchanged.It also shows that this model can optimize the layout of the enterprise according to the location safety risk to achieve the lowest overall risk value.(2) Change in the enterprise association safety risk According to Figure 5c,e, Enterprises 1 and 20 in Models 2 and 3 are arranged in the same building.To display the change in enterprise layout, the association safety risk value A comparison of Figures 5a and 11 show that with a change in the location risk of Enterprise 1, the original layout of Enterprise 1 changes from the 5th floor to the 1st floor.In Model 3, Enterprise 5 was changed from the original layout from the 5th floor to the 1st floor (see Figures 5e and 12).Moreover, the contribution of Enterprise 1 to the overall risk value of the park remains unchanged.It also shows that this model can optimize the layout of the enterprise according to the location safety risk to achieve the lowest overall risk value.
(2) Change in the enterprise association safety risk According to Figure 5c,e, Enterprises 1 and 20 in Models 2 and 3 are arranged in the same building.To display the change in enterprise layout, the association safety risk value A 1,20 between Enterprises 1 and 20 is increased from 1 to 70, and the enterprise layout changes.
As shown in Figure 13, Enterprise 1 is in Building 2, and Enterprise 20 is in Building 4 after the association safety risk is changed.After raising the associated risks between Enterprises 1 and 20, the two are distributed in different buildings after re-layout (see Figures 13 and 14).The minimum risk values of the park calculated using Models 2 and 3 are 1091 and 1285, respectively, which have changed.The total risk value calculated using Model 2 increased by 60 compared with the original, which is smaller than the increase in the association risk value (i.e., A 1, 20 ).The results of Model 3 are similar.The numerical simulation results show that the two models can effectively avoid the layout of enterprises that significantly influence each other in the same building, to reduce the possibility of the domino effect of accidents in the park.
Sustainability 2023, 15, x FOR PEER REVIEW 24 of 28  , between Enterprises 1 and 20 is increased from 1 to 70, and the enterprise layout changes.
As shown in Figure 13, Enterprise 1 is in Building 2, and Enterprise 20 is in Building 4 after the association safety risk is changed.After raising the associated risks between Enterprises 1 and 20, the two are distributed in different buildings after re-layout (see Figures 13 and 14).The minimum risk values of the park calculated using Models 2 and 3 are 1091 and 1285, respectively, which have changed.The total risk value calculated using Model 2 increased by 60 compared with the original, which is smaller than the increase in the association risk value (i.e.,  , ).The results of Model 3 are similar.The numerical simulation results show that the two models can effectively avoid the layout of enterprises that significantly influence each other in the same building, to reduce the possibility of the domino effect of accidents in the park.

Sensitivity to Problem Sizes
We evaluate the solution performances for six distinct problem sizes to confirm the efficacy of the proposed approach.Some of the parameters are shown in Table 7, and the remaining parameters are consistent with the case in Section 6.1.We record the time required to find a Pareto-optimal solution.The objective function values and required CPU time are reported in Table 7.When there are fewer than 24 enterprises, our method is incredibly effective at providing solutions (see Table 7).In addition, Model 1 is always the easiest to solve.Overall, the models become harder to solve as the size of the problem grows. , between Enterprises 1 and 20 is increased from 1 to 70, and the enterprise layout changes.
As shown in Figure 13, Enterprise 1 is in Building 2, and Enterprise 20 is in Building 4 after the association safety risk is changed.After raising the associated risks between Enterprises 1 and 20, the two are distributed in different buildings after re-layout (see Figures 13 and 14).The minimum risk values of the park calculated using Models 2 and 3 are 1091 and 1285, respectively, which have changed.The total risk value calculated using Model 2 increased by 60 compared with the original, which is smaller than the increase in the association risk value (i.e.,  , ).The results of Model 3 are similar.The numerical simulation results show that the two models can effectively avoid the layout of enterprises that significantly influence each other in the same building, to reduce the possibility of the domino effect of accidents in the park.

Sensitivity to Problem Sizes
We evaluate the solution performances for six distinct problem sizes to confirm the efficacy of the proposed approach.Some of the parameters are shown in Table 7, and the remaining parameters are consistent with the case in Section 6.1.We record the time required to find a Pareto-optimal solution.The objective function values and required CPU time are reported in Table 7.When there are fewer than 24 enterprises, our method is incredibly effective at providing solutions (see Table 7).In addition, Model 1 is always the easiest to solve.Overall, the models become harder to solve as the size of the problem grows.

Sensitivity to Problem Sizes
We evaluate the solution performances for six distinct problem sizes to confirm the efficacy of the proposed approach.Some of the parameters are shown in Table 7, and the remaining parameters are consistent with the case in Section 6.1.We record the time required to find a Pareto-optimal solution.The objective function values and required CPU time are reported in Table 7.When there are fewer than 24 enterprises, our method is incredibly effective at providing solutions (see Table 7).In addition, Model 1 is always the easiest to solve.Overall, the models become harder to solve as the size of the problem grows.

Figure 1 .
Figure 1.A schematic diagram of the distribution of enterprises in the new industrial park.

Figure 1 .
Figure 1.A schematic diagram of the distribution of enterprises in the new industrial park.

Figure 2 .
Figure 2.An illustrative example of the enterprise location safety risk.

Figure 2 .
Figure 2.An illustrative example of the enterprise location safety risk.

Figure 3 .
Figure 3.The enterprise association safety risk.

Figure 3 .
Figure 3.The enterprise association safety risk.

4. 1 .
Mathematical Model 1 4.1.1.Notations Parameters: I → Set of enterprises, indexed by i ∈ I. B → Set of buildings, indexed by b ∈ B. S → Set of floors, indexed by s ∈ S. Q bs → Size of floor s in building b.F i → Area required by enterprise i.K is → Risk value of enterprise i on floor s.Z is → Expected rent of enterprise i on floor s.Decision variables: x bis = 1 Enterprise i on floor s in building b 0 otherwise 4.1.2

Figure 5 .
Figure 5.The specific enterprise layout result obtained from the three models.Figure 5.The specific enterprise layout result obtained from the three models.

Figure 6 .
Figure 6.Variations of the Pareto frontiers with different floors.Figure 6. Variations of the Pareto frontiers with different floors.

Figure 6 .
Figure 6.Variations of the Pareto frontiers with different floors.Figure 6. Variations of the Pareto frontiers with different floors.

Figure 7 .
Figure 7.The Pareto frontier with a different number of buildings.Figure 7. The Pareto frontier with a different number of buildings.

Figure 7 . 28 (Figure 8 .
Figure 7.The Pareto frontier with a different number of buildings.Figure 7. The Pareto frontier with a different number of buildings.

Figure 8 .
Figure 8. Variations of the Pareto frontiers with different building numbers.

Figure 9 .
Figure 9.The Pareto frontier with a different number of enterprises.

Figure 11 .
Figure 11.The enterprise layout obtained from Model 1 (the location safety risk was changed).

Figure 12 .
Figure 12.The enterprise layout obtained from Model 3 (the location safety risk was changed).

Figure 11 .
Figure 11.The enterprise layout obtained from Model 1 (the location safety risk was changed).

Figure 10 .
Figure 10.Variations of the Pareto frontiers with a different number of enterprises.

Figure 11 .
Figure 11.The enterprise layout obtained from Model 1 (the location safety risk was changed).

Figure 12 .
Figure 12.The enterprise layout obtained from Model 3 (the location safety risk was changed).

Figure 12 .
Figure 12.The enterprise layout obtained from Model 3 (the location safety risk was changed).

Figure 13 .
Figure 13.The enterprise layout obtained from Model 2 (the association safety risk was changed).

Figure 14 .
Figure 14.The enterprise layout obtained from Model 3 (the association safety risk was changed).

Figure 13 .
Figure 13.The enterprise layout obtained from Model 2 (the association safety risk was changed).

Figure 13 .
Figure 13.The enterprise layout obtained from Model 2 (the association safety risk was changed).

Figure 14 .
Figure 14.The enterprise layout obtained from Model 3 (the association safety risk was changed).

Figure 14 .
Figure 14.The enterprise layout obtained from Model 3 (the association safety risk was changed).

Table 1 .
A summary of the literature pertaining to the layout optimization problem.

Table 2 .
A relational table of decision variables.

Table 3 .
The basic building information of the park.

Table 4 .
The optimal Pareto solutions of Model 1 with different floors (partial).

Table 5 .
The optimal Pareto solutions of Model 2 with different floors (partial).

Table 6 .
The optimal Pareto solutions of Model 3 with different floors (partial).