COST OPTIMIZATION OF MULTIUNIT CONSTRUCTION PROJECTS USING LINEAR PROGRAMMING AND METAHEURISTIC-BASED SIMULATED ANNEALING ALGORITHM

The article presents the cost optimization model for multiunit construction projects. Multiunit projects constitute a special case of repetitive projects. They consist in the realization of many different, when it comes to size, types of residential, commercial, industrial buildings or engineering structures. Due to the specific character of construction works, actual schedules of such projects should not only take into account real costs of construction, but also be subject to specific restrictions, e.g. deadlines for the completion of units imposed by the investor. To solve the NP-hard problem of choosing the order of units’ construction there was metaheuristic algorithm of simulated annealing used. The objective function in the presented optimization model was the total value of the project cost determined on the basis of the mathematical programming model, taking into account direct and indirect costs, costs of missing deadlines and costs of work group discontinuities. In the article, an experimental analysis of the proposed method of solving the optimization task was carried out in a model that showed high efficiency in obtaining suboptimal solutions. In addition, the operation of the proposed model has been presented on a calculation example. The results obtained in it are fully satisfying.


Introduction
Problems of scheduling construction projects have currently been the subject of many studies. Due to the method of conducting works, we usually divide them into undertakings with non-repetitive character and repetitive projects. Repetitive projects are primarily characterized by the possibility of their division into parts (units), which may be, for example, work zones, sections of a certain size (e.g. length), individual floors of buildings or entire buildings or structures. Each part of a given undertaking requires the implementation of the same or a similar set of activities using the available resources, which are construction crews. Examples of repetitive projects include single-family housing estates, groups of buildings, highrise multi-story buildings, pipeline networks, roads and highways.
In order to schedule repetitive projects there is a concept of the Line of Balance used (Arditi, Tokdemir, & Suh, 2002). On its basis, the following techniques have been developed: Linear Scheduling Method (Chrzanowski & Johnston, 1986), Repetitive Project Model (Reda, 1990), Linear Scheduling Model (Harmelink & Rowings, 1998), Repetitive Scheduling Method (Harris & Ioannou, 1998). The features of the above techniques include, among others, the ability to distinguish critical activities or their parts in the project schedule, the possibility of taking into account non-linear (centered) activities in the schedule (Harmelink & Rowings, 1998;Harris & Ioannou, 1998), or ensuring the continuity of activities (Chrzanowski & Johnston, 1986;Reda, 1990). Current research on the methods of scheduling repetitive projects most often focus on the problem of optimal creation of their schedules. In the optimization models of these projects, different parameters, constraints and accepted criteria may be taken into account. Due to their nature, while optimizing their schedules, mathematical methods are currently most commonly used (e.g. linear programming (e.g. Reda, 1990;Ipsilandis, 2007;Biruk & Jaśkowski, 2017;Radziszewska-Zielina & Sroka, 2017a, 2017b, 2018, dynamic programming (e.g. Moselhi & Hassanein, 2003), neural networks (e.g. Adeli & Karim, 1997)) or metaheuristics (e.g. genetic (Hegazy, Elhakeem, & Elbeltagi, 2004;Agrama, 2014), evolutionary (e.g. Rogalska, Bożejko, & Hejducki, 2008), hybrid (e.g. Ezeldin & Soliman, 2009), simulated annealing (e.g. Chen & Shahandashti, 2007) algorithms). They are mostly used to solve optimization problems related to, among others, the duration of the project, its cost, taking into account the condition of the imposed deadline for its completion, the total time of downtime of the construction crews carrying out the project.
As part of repetitive undertakings, multiunit projects can be distinguished, which include the construction of many residential, commercial, industrial buildings or engineering structures. The most common assumption in this type of projects is the constant and unchanging sequence of the implementation of units. This assumption may be dictated by existing organizational constraints or limitations imposed by the investor in connection with, for example, a specific sales order of the project units. However, in the general case, the order can be arbitrary. Assuming any size of works that make up a particular object, we get the opportunity to create an optimal schedule for a given multiunit construction project. In this case, the contemporary achievements of the job scheduling theory in the field of flow shop problems will be helpful (Gupta & Stafford, 2006), which are NP-hard already for n = 3 objects of the project. These problems are mainly used in electronics, the chemical industry, the automotive industry, and less frequently in construction (Bożejko, Hejducki, & Wodecki, 2012;Bożejko, Hejducki, Uchroński, & Wodecki, 2014;Podolski, 2016Podolski, , 2017. The article (Bożejko et al., 2012) presents the problem of a multiunit construction projects, in which a solution was sought to minimize costs, including only penalties for exceeding the imposed deadlines for the construction of units. The solution, which was the order of units' realization, was obtained using the scatter search algorithm. Bożejko et al. (2014) presents the problem of a multiunit construction project, in which there was the possibility of overlap of works performed in one unit, i.e. there was, in certain periods of time, the possibility of simultaneously carrying out two different works in one unit. The solution, which was also the order of units' realization, was obtained using the construction and tabu search algorithms. In the article (Podolski, 2016) a discrete problem of optimization of a two-criteria cost/ time was presented, in which the decision variables were: the order (permutation) of realization of units (buildings) and a discrete set of ways (subcontractors' offers) of works realization. Each variant included exact cost and time of the works in a given unit. The optimization task consisted in solving the single-criterion task of discrete optimization -minimizing the cost of the project with the specified limitation for the time of its realization (the sum of the selected subcontractors' offers). This problem was solved using a developed by author, two-stage, suboptimal algorithm using the simulated annealing algorithm. In the article (Podolski, 2017) a model of scheduling a multiunit project was presented, in which there is a possibility of applying more than one working group to perform a given work in the unit and order dependencies between works in a given unit in the form of a network as in the CPM method. The task of the minimization of the project duration was solved using the program developed by author in the Mathematica system, using the tabu search algorithm similarly to the flow shop problems with the parallel machines (hybrid flow shop).
Due to the specificity of construction works carried out by construction crews, actual schedules of multiunit construction projects should not only take into account the real costs of works implementation, but also be subject to specific constraints. These include conditions regarding the dates of completion of the facilities imposed by the investor, additional costs caused by the failure to maintain continuity of works performed by construction crews. The work will present a new flow shop model of a multiunit construction project, in which the actual (direct and indirect) costs of the project's implementation are taken into account, with reference to: costs of missing deadlines for construction of units, costs associated with failure to keep the continuity of works being done and variable sequence of units realization. Consideration of the aforementioned parameters and limitations required the use of a metaheuristic simulated annealing algorithm. While calculating the value of the objective function in this algorithm, which was the total cost of a multiunit construction project, the cost optimization of the schedule was made. A novelty in this article in relation to the previous ones is the modeling of cost / time dependence as a continuous function, and the assume of the actual costs of a multiunit project (direct costs of the works, indirect costs, costs of failure to meet the deadlines for the realization of works in the units, and costs of downtime for working groups) in combination with a discrete function of the order of realized units. In this article, due to the lack of research in the literature regarding the presented form of the simulated annealing algorithm, an experimental analysis of the developed metaheuristic confirming the high efficiency of the optimization algorithm was presented. The presented optimization model of a multiunit project was illustrated by a computation example.

Description of the multiunit construction project considered in the work
The multiunit construction project model considered in the work can be regarded as deterministic, when all parameters and constraints of the model are precisely defined and known from the beginning. Each of the project units requires the same set of works to be done, which are carried out in the same order. For these units, the following works will be performed successively in line with the same technology, for instance earthworks, foundation works, walls with ceilings, roof trusses. The fixed sequence of works is characteristic of buildings with a simple construction such as: single-family houses, simple industrial or engineering structures. Works in the units will be carried out by construction crews specialized to perform only one type of work, moving from the previous unit to the next one. It is assumed that the size of units included in the project will be arbitrary. They are not characterized by proportionality between the size of units and the same labor consumption of one kind of works. It follows from the above conclusion that the shape of the schedule of a multiunit project defined in this way will be influenced by the order in which the units are constructed. In the considered project model, each of the works of a given unit will be described by its normal duration and its normal cost. It is assumed that it is possible to shorten the normal time of a given work to the border time using more efficient employees, technologies, the involvement of additional workers or machines. It is connected with raising the cost of its implementation to the boundary cost. It is assumed that the dependence of the cost per time for a given work in the project site is linear. For the contractor of the project, the criterion of choosing the optimal schedule will be the total cost of the project. It will include a direct cost understood as the sum of costs of all works in project, indirect costs depending on the duration of the project, penalties for failure to comply with the deadlines for the implementation of units imposed by the investor and additional costs of the contractor due to discontinuity of works of one type in units by individual construction crews.

Optimization model of the multiunit construction project under consideration
The following optimization model of the multiunit construction project considered in the work was defined: Parameters: -The project constitutes a set of units (buildings, structures) Z = {Z 1 , Z 2 , Z 3 , ..., Z i , ..., Z n }. -In order to perform the works there are construction crews doing work of one kind and forming a set of The normal duration of the work O ij performed by the crew is tn i,j > 0. The normal duration of work O ij corresponds to the normal cost of its implementation equal to kn i,j > 0. Due to the change in the amount of resources used or the change in the technology of works, it is possible to shorten each normal time tn i,j to the value of the boundary time tb i,j > 0, where tb i,j < tn i,j , which involves raising the cost to boundary cost equaling kb i,j > kn i,j . -Indirect unit cost per day is assumed at the level of ki.
-The implementation dates of each of the n units imposed by the investor are assumed and equals Td i .
-Unit costs are assumed for failing to meet the deadlines of the individual units per day and equals kl i . -Unit costs of discontinuities (downtimes) of crews are assumed and equal kc j .
Constraints: -The order in which the works are performed resulting from the technology is assumed: -It is assumed that at any time, each crew B j can perform only one work at a time and that only one work O ij can be performed at any time in the unit Z i . -It is assumed that the O ij ∈ O i work is carried out continuously by the crew B j . For such adopted parameters and constraints there is a mathematical programming model formulated, in the following form: -Objective function: (2) Tf 1,1 ≥ t 1,1 ; (7) The objective function (1) is the sum of direct costs K d , indirect costs K i , costs related to failure to meet the deadlines for the units K l and the costs of downtime for crews K c . This cost should be minimized. Direct costs (2) are the sum of the costs of all work performed by all crews in all units depending on the duration of the work. The function of works costs depending on the duration of this activity kd i,j (t i,j ) is a decreasing linear function. In fact, the dependence of direct costs on time is a discrete dependence determined on the basis of negotiations between the subcontractor and the general contractor. The linear dependence between the time of work and its cost is a fairly good approximation. The equation of the linear function of the direct cost of a given work to its duration can be determined from the Eqn (14), while the value of the angular coefficient and the free expression from the Eqns (15) and (16).
Indirect costs (3) depend on the amount of indirect unit costs and the duration of the entire project Tf n,m . Costs related to non-compliance with the directive deadlines (4) are calculated as the sum of unit cost products for non-compliance with the deadline and the duration of this delay for all units. The unit cost for non-compliance with the directive deadline may vary for each of the units. Penalties are specified in the contract between the investor and the general contractor and usually amount to between 0.05% and 0.2% of the gross contract value for each day of delay. Costs of downtime for crews (5) are calculated as a sum of product of unit cost for downtime in crews work and duration of downtime for all crews. Duration of downtime is calculated on the basis of completion date of crew work in all units, date of termination of work on the first unit and duration of crew work on remaining units.
The decision variables in the model are: t i,j, denotes the assumed duration of work O i,j , Tf i,j denotes the deadline for completion of work O i,j , l i denotes the time of exceeding the directive deadline for individual units, and π denotes the order of executing units from the set Z: π = (π(1), π(2), ..., π(i), ..., π(n)).
The model includes the following constraints. Eqn (6) limits the duration of activities to the interval between the limit (minimum) and normal (maximum) time. Eqns (7), (8), (9) specify the assumptions regarding the order of works execution. The deadline for completion of work on the first site by the first working group must be greater than the duration of this work (7). The deadline for completing the work of the j-th crew on the i-th unit must be larger than the completion date of the work of the same crew in the previous unit enlarged by the time of duration of the crew's work in the i-th unit (8). The deadline for completing the work of the j-th crew on the i-th object must be larger than the completion date of the work of the previous crew in the same unit enlarged by the time of duration of work in this unit by the j-th crew (9). The condition (10) allows its users to determine the value of the variable l i , which determines the time of exceeding the directive deadlines. Conditions (11), (12) and (13) limit the values of decision variables to non-negative numbers.
The duration of activities and the corresponding costs of performing activities are linear in nature. For the set order of unit execution, the presented optimization model becomes a linear programming model that can be solved using, for example, the Simplex method. However, con-sidering any order π in which units are executed causes the fact that the presented model becomes an NP-hard problem of discrete optimization, which has not been examined in field of scheduling of construction project yet. The NP-hardness of the problem in the model under consideration results from the following assumption: tn i,j = tb i,j , kn i,j = kb i,j = 0, ki = 0, kl i = 0, kc j = 0, and the objective function is Tf n,m , whereas we are looking for its minimum value. Taking into consideration the above assumptions we obtain a classic permutational flow shop problem FP ||C max , which is strongly NP-hard (Gupta & Stafford, 2006). In order to find the optimal solution, the considered model requires the use of an approximate algorithm that will allow us to obtain solutions close to optimal. In this work, a metaheuristic simulated annealing algorithm was used to obtain solutions in the model.

Method for solving the optimization problem
The problem presented above in the optimization model is an NP-hard problem of discrete optimization. An important difficulty in solving it is the fact that two different groups of decision variables can be distinguished in it. The first of them is related to the selection of optimal durations of the works t i,j and the corresponding costs of performing the activities in the project taking into account the accepted objective function of the total cost. Due to the continuous and linear nature of these variables, this problem will be solved using the Simplex method, which in the considered model allows us to obtain an optimal solution. Such a solution in the form of optimal durations of activities t i,j is possible only for the determined order of executing construction units π. The second group of decision variables is represented by the order in which units are executed. The nature of this decision variable is discrete, and the problem of searching for the optimal order of performing objects is NP-hard. To solve this problem, an approximate, metaheuristic simulated annealing (SA) algorithm was chosen because of its confirmed high effectiveness in solving of optimization tasks in multiunit construction projects (Podolski, 2008). SA algorithm has been proposed in the work of Kirkpatrick, Gelatt, and Vecchi (1983). This algorithm uses analogous to the thermodynamic process of cooling the solid in order to introduce the trajectory of the search of the local extremum. States of solid matter are seen analogously as individual solution to the problem, whereas the energy of the body as the value of the objective function. During the physical process of cooling the temperature is reduced slowly in order to maintain energy balance. The SA algorithm starts with the initial solution, usually chosen at random. Then, in each iteration, according to established rules or randomly, there is solution π' selected from the base neighborhood π. It becomes the base solution in the next iteration, if the value of the objective function is better than the current base solution or if it otherwise may become the base solution with the probability of: p = exp(-D/T i ), where D = c(π') -c(π), T i -the temperature of the current iteration i, c -the objective function. In each iteration there are m draws from the neighbourhood of the current basic solution performed. The parameter called the temperature decreases in the same way as in the natural process of annealing. The most frequently adopted patterns of cooling are: N -1, T 0 -the initial temperature, T N -the final temperature, N -number of iterations, λ i -parameter. In the algorithm there are usually parameter values T 0 , T N , N adopted and parameter λ i is calculated. The relationship T 0 > T N should take place, whereas T N should be small, close to zero. Below, there is presented a general method of SA algorithm used to solve the flow shop problem.
Step 2. Change the temperature T according to a defined pattern of cooling.
Step 3. If T > T N , return to step 1, otherwise STOP.
SA algorithms are used to solve many optimization problems, including flow shop problems considered in the context of discrete optimization problems (e.g. Ogbu & Smith, 1995;Ishibuchi, Misaki, & Tanaka, 1995). Good results obtained in applications allow us to treat SA algorithms as one of the strongest tools. The optimization problem considered in the work will be solved using the algorithm presented in Figure 1. Its steps are consistent with the above-mentioned SA algorithm. When calculating the assumed objective function for a given π sequence, the problem of minimizing this function is optimally solved using the Simplex method.
In the SA algorithm used in the work, the following assumptions regarding its parameters were adopted: -the neighbourhood N(π) contains permutations generated from π using the "insert" move type, -the Boltzmann acceptance function was used, -a geometric diagram of cooling was adopted, i.e.
T i+1 = λT i , the initial temperature T 0 = 30, λ = 0.99, the number of solutions considered at a fixed temperature -300, minimum temperature T N = 0.05. The presented algorithm has been implemented in the Python programming language using the PyMathProg linear optimization package, which uses the GLPK (GNU Linear Programming Kit) engine. Figure 1. Block diagram of the applied SA algorithm Parameters: normal limes, boundary times, normal costs, boundary costs, unit indirect cost, deadlines for objects, the unit costs of failing to meet deadlines, the unit costs for works discontinuity Start Select initial permutation pi 0 Step 0: Step 1: Calculate K (π) Step 1a: assign: Step 1b: Generate a linear programming model Step 1c: Solve the model using SIMPLEX Step 1d: Step 2: Execute x times Assign k = k + l, generate π = N(π ) k-1 Step 2a: assign: Step 2b: Calculate K (π) as in step 1 min If K (π) < K (π ), then assign π = π min min SA SA Step 2c: assign: If K (π) < K (π ) or min min k-1 Step 2d: rnd() < exp(( K (π ) -K (p))/T min k-1 min then assing π = π k otherwise assign π = π k k-1 Step 3: Calculate new temperature T Step 4:

Verification of the results obtained with the use of the applied SA algorithm
For the approximate SA algorithm presented above, the authors verified the results obtained with its help. For this purpose, two groups of examples of objects were generated randomly: n = 7 objects, m = 5 works and n = 10 objects, m = 10 works. These sizes are similar to the sizes of problems that can be encountered in the practice of multiunit construction projects. There were five examples generated in each group. In each of the examples, a minimum value of the cost of a given project was sought in accordance with the model presented in the work. Each of the examples was resolved five times. The results obtained with the help of the created software are presented in Table 1. Then, the two groups of examples analyzed here were solved optimally by means of the exhaustive search algorithm (ES algorithm). The obtained results were compared with each other by calculating the average relative error PRD(SA) of the SA algorithm: where: K SA -the value of the adopted objective function obtained by means of the SA algorithm, K ES -the value of the assumed objective function obtained by means of an exhaustive search algorithm. The average relative errors of the applied SA algorithm are presented in Table 1. The values of these errors are small and do not exceed 0.5%, which confirms the high effectiveness of the SA algorithm in searching for optimal solutions in the subject.

Case study
The contractor, at the request of the investor, is to carry out a project consisting in the construction of n = 12 residential buildings (in the model called units). Each of them requires the execution of m = 9 works executed in a fixed order. The undertaking will be implemented entirely by subcontractors (in the model called construction crews). For each type of works, the contractor received from subcontractors the range of deadlines in which the work in a given building can be performed, and the corresponding costs. This is shown in Table 2. The scope of deadlines is defined by the normal tn i,j and boundary tb i,j date of the execution of a given type of work in individual buildings expressed in working days and corresponding to the normal term -normal cost kn i,j and the boundary date -boundary cost kb i,j . Each work can be executed at any time between the limit time and normal time at proportional costs. The contractor also incurs indirect costs related to running and supervising the construction in the amount of 300 € for each day of the project implementation. The contractor has a deadline for the implementation of individual buildings imposed by the investor. The directive deadlines are for subsequent objects (in working days): 160,200,240,280,320,360,400,440,480,520,560,600. For failure to meet the deadline, a penalty of 200 € is imposed for the general contractor for each day of delay. In addition, in connection with the concluded contracts with subcontractors, the contractor undertook to ensure the continuity of subcontractors' work j = 4 (roof truss, cover), j = 5 (installations), j = 7 (plaster, floors, attic), j = 9 (tiling, painting, sanitary assembly). The penalty for each day of discontinuity of work of contractors has been determined with subcontractors and amounts to 300 € for contractors 4 and 5 and 200 € for contractors 7 and 9. Additional technological and organizational limitation is that the work cannot start if the work of the same type in the previous building did not end and if the previous work in the same building was not completed. The number of possible schedules due to the order in which the buildings are realized in the example shown is 12! ≈ 4.8 •10 8 . The number of possible schedules due to the continuity of the cost-time dependency is infinite. Due to the possibility of choosing the time of works' execution, taking into account indirect costs, imposed dates of completion of individual buildings and costs to be incurred due to the lack of continuity of subcontracting work, the contractor will seek to set a schedule to minimize the total project costs. For the initial schedule of π 0 = (1,2,3,4,5,6,7,8,9,10,11,12) and assumed times of normal execution of all works, total cost of 1292.91 thousand of euro was obtained. The duration of the entire project is at the level of 625 days. The schedule for the initial π 0 rank without optimization is shown in Figure 2(a).
For the initial schedule of π 0 after solving the linear optimization task using the Simplex method and the created software, total costs of 1065.70 thousand Euro were obtained with a realization time of 457 days. This means a 17.6% reduction in costs and an acceleration of the project implementation time by 26.9% compared to the solution without the use of optimization. The schedule for the initial order π 0 after linear optimization is shown in Figure 2(b).
The next step in the search for an optimal solution was the use of the SA algorithm in order to find the optimal schedule of the project that minimizes the assumed cost goal function, taking into account the possible change in the order of construction objects execution. The stop condition for the algorithm's operation was to achieve the assumed minimum temperature T min = 0.05. The algorithm performed 191100 iterations in the SA loop, which is less than 0.04% of all possible scheduling. The run of SA algorithm is presented in Figure 3. The smallest cost of the project is 1045.28 thousand Euro obtained for the following scheduling of objects: π SA = (6,7,10,2,3,9,1,5,11,12,4,8). The deadline for the project implementation was 438 days. In relation to the optimal solution for the initial scheduling π 0 , the total cost was improved by 1.9% (20.42 thousand Euro) and the deadline for implementation was improved by 4.2% (19 days). In relation to the Figure 2. Schedules for the implementation of project in the case study: (a) -for π 0 = (1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12) without optimization; (b) -for π 0 with optimization; (c) -for π SA = (6, 7, 10, 2, 3, 9, 1, 5, 11, 12, 4, 8) with optimization using SA algorithm initial schedule (without optimization) for the schedule π 0 , the project costs were reduced by 19% and the implementation time decreased by 30%. The obtained schedule is shown in Figure 2(c).

Conclusions
Scheduling of multiunit construction projects, which is a case of repetitive projects, is a problem in which optimization tasks often occur. Due to the specific character of multiunit construction projects, decision variables in such undertakings may be continuous or discrete. An important difficulty in finding solutions introduces simultaneous occurrence of a continuous and discrete decision variables in the project optimization model, which took place in the model presented in the work. The discreteness of the decision variable -the order in which objects are realized in a multiunit project brought the presented model to the NP-hard optimization problem, i.e. the permutational flow shop problem. The search for solutions in it required the use of a metaheuristic simulated annealing algorithm with simultaneous solving the linear optimization task using the Simplex method while calculating the value of the adopted objective function -the total cost of the project. Due to the lack of research in the literature regarding the applied form of the simulated annealing algorithm, the article presents verification of the results provided with its use. It showed its high effectiveness in searching for optimal solutions in the model in question. The presented model of multiunit construction project scheduling can be used when determining the optimal work schedule for construction companies undertaking such projects. It enables to take into account the situation when companies use their own working crews to implement the works or intend to select specialized subcontractors who best meet the imposed conditions. The developed model can be easily implemented in scheduling programs and can support the planner when designing multiunit construction projects, helping to minimize the cost of the planned investment realization. Due to the development of dedicated software for solving optimization tasks, it is possible to extend the current model with additional technological and organizational constraints, additional factors affecting the cost of the project and taking into account other objective functions.

Author contributions
Michał Podolski was responsible for literature review, description of the model, partially for optimization model, description of applied SA algorithm, concept of the algorithm verification and wrote the introduction, sections 1, 3 and 4, partially sections 2 and 5, conclusions of the article. Bartłomiej Sroka was responsible for concept, computer implementation, verification and case study of the presented model, partially for optimization model, wrote partially sections 2 and 5.

Disclosure statement
Authors have not any competing financial, professional, or personal interests from other parties.