Profit optimization for multi-mode repetitive construction project with cash flows using metaheuristics

The article presents the profit optimization model for multi-unit construction projects. Such projects constitute a special case of repetitive projects and are common in residential, commercial, and industrial construction projects. Due to the specific character of construction works, schedules of such projects should take into account many different aspects, including durations and costs of construction works, the possibility of selecting alternative execution modes, and 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 and the best variants of works, the authors used metaheuristic algorithms (simulated annealing and genetic search). The objective function in the presented optimization model was the total profit of the contractor determined on the basis of the mathematical programming model. This model takes into account monthly cash flows subject to direct and indirect costs, penalties for missing deadlines, costs of work group discontinuities, and borrowing losses. The presented problem is very important for maintaining a good financial condition of the enterprise carrying out construction projects. 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
Scheduling is crucial for construction project planning. Properly prepared schedules allow construction managers to prioritize and allocate time for different tasks; represent the dependencies between them; manage resources, both renewable (labor, machinery, materials) and non-renewable (funds). To create a schedule, one must determine the order of tasks and their start times. This must be performed in such a way that specific objectives of the project are achieved whilst at the same time, organizational conditions are fulfilled, and planning constraints are applied. Some of the construction project planning objectives include minimizing the duration of the project; minimizing or providing evenness of the resources' consumption; maximizing the construction project economic value for the contractor. Organizational conditions may include, for example, the need to ensure the continuity of work performed by various groups of workers, or the possibility of using different construction activities performance modes. Basic planning constraints are linked to the availability of the time, and resources needed to implement them and technological dependencies between the activities.
To meet the needs of the construction market, as well as to ensure flexibility in decision-making, the authors decided to combine concepts of flow shop and multi-mode resourceconstrained project scheduling problem. The article presents a new optimization model of a multi-unit construction project with specific assumptions. Its most important assumptions are the need to perform many works in the undertaking for many different building structures, the availability of many different contractors for works in the facilities. The new scheduling model for a multi-unit construction project is illustrated by a computational example. The new presented model takes into account discounted cash flows, the knowledge of which is necessary for construction companies to make strategic decisions, and to survive on the market. In the model presented in the article, the whole project timespan is divided into n-day intervals representing billing periods (usually a month). Direct and indirect costs are calculated proportionally in each period and are discounted every n-days. Income was modeled in a form of Progress Payments (PP) which is common in the construction industry. The client pays contractors based on a monthly invoice, and report on completed works. The payments correspond to the direct and indirect costs plus the contractor's profit (and are also discounted). However, taking into account specifics of the construction industry, the authors assumed that the payments are made with a 1-month lag. To reflect the market characteristics, contractual penalties for exceeding the deadline were introduced into the model. Penalties related to the discontinuity of work of the work teams were also applied. These penalties are not discounted and are accounted for with a lag equal to one billing period. Accumulated cash flows are calculated for each n-day period. If the accumulated cash flow in a given period is negative, then we charge an additional penalty that corresponds to the cost of the loan required for the financing of the contractor's operations. The objective function assumes the maximization of the total profit-cumulated cash flow on the last day of the project's implementation (one period after completion of construction works).

Literature review
Theoretical studies of project scheduling problems are currently focused on searching for optimal schedules considering existing constraints. These are usually NP-hard optimization problems, as presented in [2,3]. Due to the variety of possible constraints and schedules' objective functions, these problems can be grouped into different categories-models of project scheduling problems (PSP). Deterministic approach for solving PSP optimization is the most common one, due to its practical aspect. The state-ofthe-art review of PSP deterministic models was presented in, e.g., [4][5][6]. One of the most common PSP models investigated by researchers is the well-known Resource-Constrained Project Scheduling Problem (RCPSP), e.g., [7]. This phenomenon results from the wide possibilities of applying this problem in the practice of industrial production.
Nowadays, in practice, more and more contracts focus on delivering non-standard products, or ones that were individually agreed on with the client. That is why scheduling production processes are often in line with the principles of project management. The same principle is also characteristic of the construction industry. A generalization of RCPSP is the Multi-Mode RCPSP (MRCPSP or rarely MMRCPSP). In such problems, each activity can be executed in one of the several modes. Each mode has a specific duration and specific resource requirements [8]. Due to the introduction of modes, decision makers are able to study different projects' variants; for example, testing how allocating additional resources to different tasks will affect project duration. It is worth mentioning that with the introduction of the additional decision variables (activities' modes), the possible solutions' space (and, at the same time, the amount of time required to solve the problem) increases. In other words, the computational time required for solving a MRCPSP is longer than that of a similar RCPSP without multiple modes [9].
Generalized RCPS/MRCPS problem is obtained by replacing the makespan minimization with other agentsany regular measure of performance: different types of tradeoffs, objective functions, constraints, and conditions, e.g., total cost [10], NPV [11], quality [12], deviations from average employment level [13], total project delay [14], monthly cash demand [15], cost minimization in regard to the baseplan [16], and schedule robustness [17]. Such problems can be referred to as the Generalized Resource-Constrained Project Scheduling Problem (GRCPSP) [9]. Such problems can be divided into different categories, for example, the dependencies between time and cost of a project are taken into account in Multi-Mode Resource Constrained, Discrete Time-Cost Trade-Off problems (MRC-DTCTP) [13] which in construction are often called simply Time-Cost Trade-Off (TCT) problems. Other variations of the problem include P-MRCPSP (P at the beginning stands for Pre-emptive) in which activities can be pre-empted at any point in time and restarted at no additional cost [8] or MRCPSP with Discounted Cash Flows (MRCPSPDCF) [11], which focuses on maximization of NPV.
The procedures used to solve MRCPSP (and scheduling problems in general) can be classified as: exact, heuristic, and metaheuristic [1,3]. The exact procedures include, among others, linear programming (LP), dynamic programming (DP), and Branch and Bound method (B&B). The heuristic methods include priority rule-based heuristics [9].
Practical scheduling problems in construction can be easily qualified as NP-hard (non-deterministic, polynomial-time hard) problems. The time needed for solving such problems grows exponentially with the increase of the problem's size [2,3,9]-therefore, mathematical and heuristic methods often do not enable finding solutions to complicated construction problems within an acceptable period of time. In the view of many authors, metaheuristic algorithms seem to be the most appropriate measures for scheduling and task sequencing [1,[8][9][10][11]18].
The metaheuristic approach does not guarantee finding the optimal solution and the obtained results are often subject to their input parameters; however, they seem perfect for solving complicated, NP-hard class problems because they enable computing suboptimal (acceptable) solutions within an acceptable time frame.
A great variety of metaheuristic algorithms can be used for solving GRCPSPs, e.g., Genetic Algorithms (GA) [8,10], Simulated Annealing (SA) [11], Tabu Search (TS) [11,19], Particle Swarm Optimization (PSO) [20], Ant Colony Optimization (ACO) [21] or hybrid algorithms [22,23]. However, only a handful of conducted research activities, concerning multi-mode problems, revolve around the construction industry. Some of these cases involve real-life case studies. Özdamar and Dündar [24] optimized NPV for apartment building construction, with different modes' duration (however, activities' cost parameters were fixed). Chen and Weng [25], and Ghoddousi et al. [13] used GA to optimize time-cost trade-off for a simplified warehouse construction project. Zhang and Xu [26] minimized makespan (at the same time maximizing quality) of the hydropower plant using PSO algorithm. The same case has been optimized in terms of time, cost and quality, by Xu and Feng, with the use of hybrid algorithm [27]. Kulejewski and Rosłon reinforced the tabu search algorithm with artificial neural networks (ANN) to minimize the maximum monthly cash demand of apartment building construction [15]. Hegazy with other contributors minimized the total cost of multi-site [28] and repetitive projects [29]. Zhang et al. minimized duration for an example bridge project, their method was developed to facilitate repetitive project scheduling [30].
Methods of repetitive projects scheduling play important role in the current research regarding the construction project scheduling. Such projects can be easily partitioned into different units which can take a form of work zones, whole storys, building structures, or sections of construction objects characterized by the length (for example pipelines). Examples of the projects realized in a flow organization system include multi-story buildings, housing estates or groups of buildings, pipelines, and roads. Each unit of the project requires specialized work crews to perform assumed tasks. These work crews are moving from one unit to another [31]. A special case of repetitive projects is multi-unit projects, which include the implementation of residential, service, industrial or engineering constructions. In general, a feature of this type of project is the ability to set any order for the implementation of the works. Assuming any size of works that make up the implementation of a specific object in a multi-unit project, we get the opportunity to create the optimal schedule for the project. For n = 3 objects, optimization tasks in this problem are NP-hard. To describe and solve such problems, we can use the flow shop theory of scheduling [32]. Research on the problem of multi-unit projects scheduling is focusing mainly on the improvement of current optimization models and improvement of optimization methods [33][34][35][36][37]. In [33], the problem of scheduling with minimization of penalties for exceeding the project's objects completion deadlines is considered and solved by the metaheuristic scatter search algorithm. In [34], the multiunit project scheduling model takes into account the possibility of overlapping works in the facilities and is solved using the tabu search algorithm. The optimization criterion is the minimization of the project duration. In the article [35], the NP-difficult problem of scheduling a construction project was considered with the criterion of the sum of penalties for exceeding the deadline for building construction. The parameters of this project were represented by fuzzy numbers or random variables with a normal distribution or the Erlang distribution. In [36], a scheduling model was presented with the possibility of performing one type of work by more than one working group and with sequence relationships given by any graph. The optimization task in this model was solved using the tabu search algorithm. In [37], the scheduling model for a multi-unit project assumes a linear relationship between the time and the cost of carrying out the activity. The criterion of minimizing 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 workgroup discontinuities. The optimization task in this model was solved using the modified simulated annealing algorithm.

Model of multi-unit project
The paper deals with the issue of construction project's profit maximization. Using the terminology of the scheduling theory, this problem can be described as follows: a set of n indivisible tasks is given (in the case of this paper-construction objects) that must be performed with m machines (in this paper-teams of working groups). Each of the n construction objects requires m types of operations (in this paper-works) to be performed by m teams of working groups. Out of m teams of workgroups, only one workgroup must be selected to perform a given operation (work) in a given construction object.
It is assumed that the order of execution of construction objects for each of the teams of working groups is the same, i.e., it is a permutational problem. We also assume that each working group from among working group teams has a strictly defined duration t and cost c of the work. Between works on a given construction object, there may be technological breaks between works or there may be a partial overlap of works. Planning a schedule for carrying out such an undertaking, it is crucial to find its optimal schedule that takes into account the criterion adopted by the schedule planner (decision maker).
Finding the schedule of the project makes it possible to determine the order of execution of the objects and the set of choices for the methods of execution of works (the selection of working groups from teams of working groups). The criterion adopted in the project is the maximization of the total profit achieved in the project, more precisely the cumulative cash flow on the last day of the project implementation (an example of the cumulative cash flow chart, prepared in accordance with the principles set out in [37], is shown in Fig. 1-it takes into account the costs incurred during project's implementation: direct costs of working groups, indirect costs, contractual penalties for exceeding the deadline for completing the construction, penalties related to breaks in the work of working groups, and loan costs).
The presented problem is a generalization of the classical permutational flow shop problem, which is shown schematically in Fig. 2. In the terminology of the scheduling theory, it is denoted as the F||C max problem, according to Graham's notation, presented in [39]. In Fig. 2, the variable S denotes the input sequence (permutation) of n tasks to the system of m machines. The cardinality of the set of possible solutions is n!.

Optimization model
The optimization model of the above described problem is as follows: Parameters: • The project consists of a set of building units • To carry out the project works the teams of working groups perform one job of one type which constitute the set The cost of the work c ijk is determined by calculation of the cost of execution of the work O ij by the working group B jk located in the resources of the contractor. It may also be the cost offer of the execution of the work O ij made by the subcontractor represented by the working group B jk . It is assumed that the time of execution of works t ijk are convex, decreasing cost functions c ijk . • There is the possibility of technological gaps between the works and the simultaneous operation of multiple working groups in the units assumed. Durations of intervals between a given work and the next work (s ij ≥ 0) or the length of the simultaneous duration of a given work and the next work (s ij < 0) in the unit for a set of works O i are given in vector s i = [s i1 , s i2 , s i3 ,…,s ik ,…,s im ]. These times should be understood as the minimum constraint and can take any value. In the further work, there will be called couplings between units. Constraints: • The order of execution of the works resulting from work technology is assumed such that • It is assumed that each working group from the team B j can perform only one job at a time. • It is assumed that the work O ij ∈ O i is performed continuously by the working group B jk ⊂ B j in time t ijk > 0.
Decision variables are the order π of execution of units, which for each of the working group, is the same and takes the form of a permutation and a set of numbers of ways of work execution (from k = 1 to k = p) in all units of the project is , R ij is the number of ways of realization (from k = 1 to k = p) of the work j in unit i, and R is the set of all possible ways to carry out the works in the project. The value of the number k of the way of realization R ij enables allocation of the working group B jk from the team B j to the work j in the object i. The form of decision variable R uniquely identifies the allocation of working groups to realization of works in the units. Therefore, using the decision variable R, there are uniquely established not only the durations of individual works carried out in the units but also their cost. After the adoption of the decision variable, R the durations of works t j from the set O i is as follows: where t ij is the duration of the execution of the work j in the object (unit) i. Similarly, in consequence of the adoption of the decision variable R, the set of costs c i of works from the set O i is as follows: where c ij is the cost of implementing work j in the unit i.
The deadlines for the individual works for the decision variables π and R can be determined from the recursive formula: The duration of the entire project F n,m (time execution of all works in the units) for π* ∈ П and for the decision variable R * ∈ R is F n,m (π * , R * ) = F π*(n),m.
The deadlines for the performance of individual works and their cost can be found in time O(nm). The number of possible solutions to the presented model is n!×p mn .
Objective function will be cumulated cash flow in the last billing period of the project which represents the total profit for the contractor. The cash flow will be maximized. Variables in the objective function are.
• F i,j is the finish time for works performed by a work crew j on unit i, • F n,m is the finish time for construction works, • h is the billing period, h ∈ {1, … , H} , where H is the last billing period, • TI is the time interval (in days), which represents the duration of each billing period, • t i,j is the duration of activities performed by a work crew j to finish all required works on unit i, For such adopted parameters, constraints and variables, there is a mathematical programming model formulated, in the following form: The objective function (2) maximizes cumulated cash flow in the last billing period of the project which represents the total profit (TP) for the contractor. The number of billing periods (3) is calculated on the base of the time interval and finish time for construction works. Equation (4) presents the iterative method for cumulated cash flow calculation in the following billing periods.
Negative cash flow results in a necessity for finding additional funds for the conduction of the project's works. For contractors, this means taking interest-bearing loans or losing the option of investing their own funds (lost profits). Binary variable (5) is modeling this phenomenon by granting penalties in each billing period with negative cash flow.
(2) CF H+max (dInc;dPen) → max Production cost (6) in each period is a discounted sum of indirect (7) and direct (8) costs in this period. Production value is invoiced after completion of a billing period and is received (income) by the contractor from the client with an agreed delay ( d Inc ) . Production value is subject to discount. Penalties for a delay of works (10) and discontinuities of work crews (11) are also accounted with a delay. First types of 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. Downtime penalties are related to the contractor's need to keep the brigades ready for work at all times.
The described model can be represented in the form of a disjunctive graph G(π), an example of which, with highlighted critical path is shown in Fig. 3. The form of the graph is dependent upon an established decisive variable π: G(π) = (N, E(π)),where N is a set of nodes, E(π)-set of arcs. It is assumed that N = {1,.., i,…., n} {1,…, j,…, m} is a set of nodes representing works in units of project. Weight of node (i, j) is equal to the time of work execution t π(i),j . The set of E(π) = E F E S (π) depends on the assumed decisive variable π. Horizontal arcs (sequential, representing processing order of units) from the set E S (π) are between nodes π(i-1) and π(i), where i = 1,…, n. Vertical arcs (technological) from the set E F are between node standing for work j and j-1. The weight of the vertical arc from set E F is the coupling between units s π(i),j . The above presented model is NP-hard optimization problem, because of assumptions from permutation flow shop problem (problem F||C max ), which is strongly NP-hard. To solve the optimization task, there is individual algorithm proposed. This algorithm will use approximate metaheuristic simulated annealing algorithm or genetic search algorithm.

Optimization methods for the presented multi-unit model
The optimization model for the multi-unit project presented above is an NP-hard problem of discrete optimization. There are two different decision variables affecting the value of the objective function. This causes the need to construct an effective individual algorithm for solving discrete optimization tasks in the presented problem. The first decision variable is the order of execution of the project objects (units), which is represented by permutation of the execution of the objects π. This decision variable is discrete, and the issue of searching for the optimal order of execution of the units is NP-hard. Due to the small number of units (most often the maximum number of objects is a dozen or so) that make up the multi-unit project, it is proposed to use the simulated annealing algorithm to solve the sequencing optimization task.
The second decision variable is a set of numbers presenting modes of performing works by available contractors. This decision variable generates a much larger area of possible solutions depending on the number of objects, types of works and the number of modes of performing works for each element of a given unit. Therefore, in this article, it is proposed that this decision variable will be solved using a suitably adapted simulated annealing algorithm or genetic search algorithm.
The method of solving the discrete optimization task for the model of a multi-unit project presented in chapter 3, taking into account both different decision variables, is presented in the form of block diagrams in Figs. 4 and 5. Its steps are consistent with SA algorithm presented in chapter 4.1. When calculating the adopted objective function for a given order π, the problem

Simulated annealing algorithm
Simulated annealing (SA) algorithm has been proposed in the work of Kirkpatrick [40]. This algorithm uses analogous to the thermodynamic process of cooling the solid 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 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 where π = c(π)c(π´), T i is the temperature of the current iteration i, c is the objective function. In each iteration, there are m draws from the neighborhood 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 geometrical or logarithmic ones. Below is the general algorithm of the SA method used to solve flow problems, which will be used to solve the problem under consideration.
SA algorithms are used to solve many optimization problems, including flow shop problems which are considered in the context of discrete optimization problems.

Genetic search algorithm
Genetic search algorithms (GS) use the principles of evolution in nature, which lead to the best adaptation (optimization) of individuals to the conditions found in a given environment. GS algorithm concept was presented in [41,42]. They use a population of individuals (solutions), which are then processed during the selection, actions induced by the use of genetic operators, and the survival phase. The population in GS algorithms is a set of individuals representing solutions. Each solution is coded by a set of attributes stored in the genetic material (chromosomes, genes). There are many coding methods specific for various optimization problems, e.g., for flow shop problems, subsequent solutions are coded in the chromosome directly using permutation [43]. The population is processed by means of cyclically consecutive processes: reproduction, crossing, and mutation as well as survival or selection. In the reproduction phase, individuals are reproduced in proportion to the measure of adaptation to the environment. The adaptation function, which is a measure of adaptation, can be, e.g., the value of the objective function for a given solution. This process means that individuals with better adaptation will have more descendants in the next generations. Individuals selected from the population form the so-called parental pool from which the pairs are selected (so-called parents) providing individuals of a new generation. They are created using the genetic crossing operator. Then, the mutation process is carried out causing changes in the genetic material, which usually occur with a low probability allowing for the slow introduction of innovation in a generation. In the survival (selection) phase, individuals are selected that will form part of the new population.
Selection can be carried out, e.g., in accordance with the roulette principle (rank selection), in which better chances are given to better-adapted individuals or on a tournament basis. The conditions for stopping the algorithm can be: time limit, the maximum number of iterations, achieving a satisfactory value of the objective function, or the optimal value. The GS method contains many elements that can be freely defined: chromosome coding method, crossover and mutation operators, adaptation function, parent pool selection, and parent matching scheme, survival scheme. The following is the general algorithm for the GS method: GS algorithms are currently successfully used to solve a number of optimization problems, including those related to the theory of task scheduling [43]. The GS method has imperfections that can manifest, mainly early convergence to the local extreme or poor convergence to optimal solutions or close to them. However, it behaves quite well for examples of small sizes, such as the model presented in the article.

The SA + SA form of solving the problem according to the presented model
The method of solving the optimization task in the presented model is consistent with the SA algorithm steps. When searching the set of possible sequencing (permutations) of object execution (the first decision variable), the following assumptions were made regarding its parameters: • the N( ) environment contains permutations generated from π using the "replace" motion, • Boltzmann acceptance function was used, • the geometric cooling scheme was adopted, i.e., T i+1 = T i and the initial temperature T 0 = 60 , = 0.7 , the number of solutions considered at the set temperature: 2, the minimum temperature T min = 0.01.
While calculating the value of the objective function, for a given permutation π, the task of optimizing profit maximization is solved, taking into account the second decision variable, i.e., a set of numbers presenting ways of performing works (modes). For this stage, the following assumptions were made regarding the parameters of the SA algorithm used in it: • the neighborhood of a given set of numbers presenting the ways of performing works contains sets generated from the output set by means of the "change to another random" movement, • Boltzmann acceptance function was used, • the geometric cooling scheme was adopted, i.e., T i+1 = T i and the initial temperature T 0 = 60 , = 0.99 , the number of solutions considered at the set temperature: 6, the mum temperature T min = 0.01.
The presented approach using the SA algorithm for both decision variables will be referred to as the SA + SA algorithm in the following part of the article and is presented in the form of a block diagram in Fig. 4. The SA + SA algorithm has been implemented in the Python programming language.

The SA + GS form of solving the problem according to the presented model
Similarly to the SA + SA algorithm, the method of solving the optimization task in the presented model is consistent with the SA algorithm steps. While searching the set of possible ordering (permutations) of object execution (the first decision variable), assumptions regarding its parameters were made the same as in the SA + SA algorithm presented in chapter 4.3.
In contrast to the SA + SA algorithm, the GS algorithm was used only when calculating the value of the objective function, when for a given permutation π, the task of optimizing profit maximization is solved taking into account the second decision variable, i.e., a set of numbers presenting the ways (modes) of performing works. For the first decision variable, i.e., the order of execution of the units (objects), it is assumed that the SA algorithm will still be used. Further in the article, the algorithm name SA + GS will be used for this combination of algorithms. The following assumptions were made regarding the form and parameters of the GS algorithm in the SA + GS algorithm: • a randomly created population consists of individualssets of numbers representing methods of works implementation R, its size is equal to 90, • selection-parents are selected using tournament selection, 40% of individuals are selected, • crossing-the single point crossing operator was used, parents are chosen randomly, • mutation is implemented on individuals of the generation by means of the "change by 1" motion applied to a randomly selected individual from the generation, part of the population is subject to mutation: 5%, probability of mutation: 0.01. • the maximum number of iterations of the GS algorithm is 90.
The SA + GS algorithm has also been implemented in the Python programming language and is presented in the form of a block diagram in Fig. 5.

Verification of the results
The authors verified the results obtained with the methods presented above. For this purpose, three groups of examples of objects were generated randomly: n = 6 units with m = 2 works, n = 5 units with m = 3 works, and n = 4 units with m = 4 works. Such sizes are characteristic for small problems in the practice of multi-unit construction projects. There were five examples generated in each group. In each of the examples, a maximum value of the total profit of a given project was sought in accordance with the presented model. Each of the studied examples was resolved five times. The results obtained with the help of the original software are presented in Table 1. Then, the three groups of analyzed examples were solved optimally by means of the exhaustive search (ES) algorithm. The obtained results were compared with each other by calculating the percentage relative difference PRD(SA + SA) of the SA + SA algorithm and PRD(SA + GS) of the SA + GS algorithms: where K SA+SA is the value of the adopted objective function obtained by means of the SA + SA algorithm, K SA+GS is the value of the adopted objective function obtained by means of the SA + GS algorithm, and K ES is the value of the assumed objective function obtained by means of an exhaustive search algorithm. The average PRD of the applied algorithms are presented in Table 1. The values of these differences are quite small, which confirms the high effectiveness of the algorithms in searching for optimal solutions in the subject. However, the results obtained by the SA + GS algorithm for small examples were slightly better with a average PRD difference of 0.452.
Additional tests were conducted for the bigger examples with n = 10, 15, 20, and m = 3, 5, 7. There were five examples generated for each size. Each of the studied examples was resolved five times. The number of examples for each of the considered sizes (in total 60 examples for all sizes) and the number of attempts to solve them is sufficient to compile statistical data, and draw conclusions about the quality of the results achieved by the tested algorithms. The average results obtained with the help of the created software are presented in Table 2. These examples were too big to be calculated by ES algorithm in a reasonable time, so the results obtained by SA + SA and SA + GS algorithms were tested using random search algorithm (RS algorithm) [44]. PRD(SA + SA) of the SA + SA algorithm and PRD(SA + GS) of the SA + GS algorithms were calculated as follows: where K RS is the value of the adopted objective function obtained by means of the RS algorithm.
The average percentage relative differences of the applied algorithms for bigger examples are presented in Table 2. In the conducted experiments, for sizes, n = 10, 15, 20, and m = 3, 5, 7 the results of the SA + GS algorithm were always better (average PRD of − 8.02%) than the results of the SA + SA algorithm (average PRD of − 1.46%). That is why SA + GS algorithm was selected for the analysis of real-life problem instances. The results of the SA + GS and SA + SA algorithms were always better than the results of the RS algorithm.

Case study
The contractor, at the request of the investor, is to carry out a project consisting in the construction of n = 5 residential buildings (units). Each of them requires execution of m = 5 works carried out in a fixed order. Technological and organizational limitation ensures that a 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 project will be implemented in full by the subcontractors. The contractor received numerous bids from the subcontractors which results in a possibility of executing each type of work in one of the three different variants (modes). Each mode has a fixed execution time of a given type of work (expressed in working days) and implementation cost (expressed in EUR). Execution modes are presented in Table 3. The calculation example was created on the basis of a model residential building (determined on the basis of the price catalogue), for which the duration of works and their direct costs were determined. On the basis of these data, the authors generated: the duration of works, and the direct cost of the object, the so-called medium building. The time and direct cost of each job were modified by a random number from − 30 to 30%. The remaining objects were generated on the basis of the so-called medium building using the uniform distribution. When creating the data on the time and cost of works, the principle was followed that any reduction in the duration of an activity resulted in an increase in its cost. The example was generated using randomness due to the inability to obtain real data for such an undertaking with sizes n = 5 objects and m = 5 works. For computation purposes, the authors used: Intel Core i5-4440, 4 GB RAM, OS Windows 10, Python 3.6 (with the uniform distribution based on numpy.random.randint function).
Between the works realized in the technological order, there are couplings between units that have been established on the basis of existing technological constraints which are shown in Table 4.
The duration of a billing period is equal to 20 working days which corresponds with one working month. This duration is set in accordance with construction practice in which works are invoiced on a monthly basis. Negative monthly cash flows force the contractor to take loans at 9% per annum (paid off every month). The indirect costs incurred by the contractor (related to running and supervising the construction) are estimated to be 730 € for each day of the project implementation. The production value is subject to discount each month at the level of 8% per annum. The contractor assumed profit at the level of 12%. The contractor has a deadline for the implementation of individual buildings imposed by the investor. The directive deadlines are 150, 210, 270, 350, and 380 (in working days) for subsequent units (objects). For failure to meet deadlines, penalties are imposed for the general contractor for each day of delay. For each building, the penalties are respectively: 540, 490, 580, 590, and 570 €. The penalties for discontinuities of work crews have been determined with subcontractors and amount respectively for each type of work (j from 1 to 5) to 200, 200, 100, 0, and 300 € for each day of stoppage. Payment and penalties delays are both equal to 20 business days.
The number of possible solutions for this problem instance is 5!×3 5×5 ≈ 1.02 ×10 14 . The optimization goal for the contractor is to maximize the total profit-TP (cumulated cash flow for the last month of the project). For the initial schedule of π 0 = (1,2,3,4,5) and assumed times and costs of execution of all works in mode 2. Total profit of 143.87 thousand of euro was obtained. The duration of the entire project is at the level of 351 days. However, not all the deadlines for units were met and some downtime penalties were charged. The schedule for the initial π 0 = (1,2,3,4,5) rank without optimization is shown in Fig. 6a. The cumulated cash flow for this schedule is presented in Fig. 7a.
The next step in the search for an optimal solution was the use of the SA + GS algorithm to find the optimal schedule of the project that maximizes the assumed goal function, taking into account the possible changes in the order of construction objects execution, and selected execution modes. The algorithm performed 51 iterations according to the parameters described in Sect. 3.
The biggest TP calculated for the project is 204.58 thousand Euro obtained for the following scheduling of units: π SA+GS = (2,3,5,1,4). The selected execution modes are presented in Table 5. The deadline for the project implementation was 283 days. In relation to the initial schedule π 0 , the total profit was improved by 42% (60.71 thousand Euro) and the deadline for implementation was improved by 19% (68 days). The obtained schedule is shown in Fig. 6b. The cumulated cash flow for this schedule is presented in Fig. 7b.

Summary and conclusions
Scheduling of multi-unit construction projects, which is a case of repetitive projects, is a problem in which discrete optimization tasks often occur. These problems are usually NP-hard due to the fact that they may be qualified as the permutation flow shop problems. In the presented model Table 4 Couplings between units s occurring between the works j and j + 1 for n = 5 units  6 a Schedule for the implementation of project in the case study for π 0 . b Schedule for the implementation of project in the case study for π SA+GS with optimization using SA + GS algorithm application of the additional parameter namely, the selection of the works' execution mode resulted in a significant increase in the number of possible solutions. Consequently, the analyzed discrete optimization problem had two different decision variables. The developed model considered monthly cash flows subject to direct and indirect costs, penalties for missing deadlines, costs of work group discontinuities, and borrowing losses. The latter factor is often overlooked in research. Meanwhile, it can have a significant impact on the outcome of projects. Especially in construction projects in which payment delays and contractual penalties occur (both aspects were included in the model). Two approaches were tested to solve the problem, a combination of two SA algorithms, and SA algorithm supported by GS algorithm. The results obtained using the second combination were slightly better than results obtained using only SA algorithms. Due to the lack of research in the literature regarding the applied form of the algorithms, 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 multi-unit construction project scheduling can be used when determining the optimal work schedule for construction companies using the flow organization system of work. It can be used in 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 multi-unit construction projects, helping to maximize the profit of the planned undertaking. Fig. 7 a Cumulated cash flows for the implementation of project in the case study for π 0 . b Cumulated cash flows for the implementation of project in the case study for π SA+GS with optimization using SA + GS algorithm Table 5 Best execution modes for π SA+GS = (2, 3, 5, 1, 4) with optimization using SA + GS algorithm

Modes
Units i = 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.
Funding The paper was partially supported by the National Science Centre of Poland, Grant OPUS no. 2017/25/B/ST7/02181.

Conflict of interest
The authors declare no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.