An Imitation and Heuristic Method for Scheduling with Subcontracted Resources

: A scheduling problem with subcontracted resources is widely spread and is associated with the distribution of limited renewable and non-renewable resources, both own and subcontracted ones based on the work’s due dates and the earliest start time. Scheduling’s goal is to reduce the cost of the subcontracted resources. In the paper, application of a few scheduling methods based on scheduling theory and the optimization algorithm is considered; limitations of these methods’ application are highlighted. It is shown that the use of simulation modeling with heuristic rules for allocation of the renewable resources makes it possible to overcome the identiﬁed limitations. A new imitation and heuristic method for solving the assigned scheduling problem is proposed. The comparison of the new method with existing ones in terms of the quality of the found solution and performance of the methods is carried out. A case study is presented that allowed a four-fold reduction of the overall subcontracted resources cost in a real project portfolio.


Introduction
The problem of business processes and project works scheduling is one of the key problems of organizational system control. Organizational systems are widely spread: examples are enterprises of various industries, multiservice communication networks, and project organizations.
According to the machine environment, scheduling problems can be divided into the main following types with the type notation in brackets [15]: identical parallel machines (Pm) where work can be performed at any of the machines [1][2][3][4][5]8,9,12,14]; Job shop (Jm) where each work has to be processed on each of the machines and all works have different routes [6,7,11,13]; Open shop (Om) where each work has to be processed on each of the machines and some of this processing time may be zero [10].
In general, the scheduling problem is a problem of work's sequence definition satisfying the given restrictions and minimizing makespan. A schedule, or calendar plan, is a set of work start dates. Commonly, researchers have investigated the renewable restricted resources that are occupied while a work is being performed and released when the work is over, as in [1,2,[4][5][6][8][9][10][11][12]14]. Examples of the renewable resources are personnel, aggregates, vehicles, and other. At the same time, there is an urgent problem to consider non-renewable resources' distribution, as in [3,7,13]. The non-renewable resources are consumed at the work input and produced at the work output. Examples are raw materials, fuel, finance, and goods.
We consider the scheduling problem on parallel machines in the presence of the orders' earliest start times and due dates. We have extended the problem by using subcontracted resources' cost optimization and accounting for restricted non-renewable resources. The optimization of subcontracted resources is relevant either for enterprises with a relatively small set of their own resources, or for enterprises with resources of constrained competence. These enterprises flexibly respond to demand fluctuations in products manufactured or services provided and attract subcontracted resources when it is necessary. Examples are construction companies and project organizations. They need to reduce the engaged subcontractors' cost while keeping with existing time restrictions to reduce the company's waste.
The rest of the paper is organized as follows. In Section 2, a literature overview is presented. Notations used are presented in Section 3. In Section 4, a scheduling problem is formulated considering own and subcontracted renewable and non-renewable resources. Application of the scheduling methods to the trial problem and their shortcomings are given in Section 5. In Section 6, we propose an imitation and heuristic method for the problem considered. In Sections 7 and 8, a case study is presented providing a comparison of the application results of the methods considered. We conclude this research and propose directions for the further work in Section 9.

Literature Overview
Network methods for the scheduling problem are introduced in [1,2]. They are intended to determine a critical path and backup time for a work. A PERT method [1] is used when works have a probabilistic duration with due dates. A GERT method [2] extends the PERT method by considering the probabilities of individual works' implementation. The methods are very useful in the case of precedence relation occurrence. Network methods do not support non-renewable resources' implementation and restrictions on the works' earliest start time.
An approximate algorithm proposed in [3] is intended to solve the scheduling problem with restricted renewable and non-renewable resources in the presence of a deadline and restrictions on the earliest start time for a work. The algorithm has two stages. At the first stage, a feasible schedule is calculated based on an assumption that all resources are non-renewable. At the second stage, the renewable resources' constraints are added to the schedule found and then works are packed to satisfy all resource constraints and minimize the execution time. A disadvantage of the algorithm is that the arrangement solution is excluded from consideration if contradictions arise between the resource constraints and work deadlines. For example, excess availability of the own renewable resources for works on the critical path can exist. In this case, companies can use subcontracted renewable resources, the cost of which needs to be optimized.
For a scheduling optimization problem in a parallel system with identical nonrenewable resources and orders' earliest start dates, two algorithms of the scheduling theory are presented in [4,5]. The algorithms are discussed in detail in Section 5.
A decentralized scheduling approach applied to a job-shop scheduling problem via agent-based simulation is presented in [6,7]. Multi-agent simulation is a popular technique Mathematics 2021, 9, 2098 3 of 22 intended to represent decision makers as a community of interacting agents [16]. In this case, each renewable resource and each single work is represented by an agent. The agents interact with each other and allocate the resources to the works by mean of negotiations. In [6], a new control algorithm based on dispatching rules is described. This algorithm is used to assign a score to each agent-resource proposal to agent-work, then the agent-work chooses the highest scored one. The resource score depends on the work's minimum processing time provided by the current resource compared to all the available alternatives, and the average cycle time of the system. An advantage of the algorithm is simultaneous multi-objective optimization of the system's throughput and cycle time, with weighted linear convolution of the two into a single objective function. In [7], the schedule is built via self-organization of software agents forming two networks of demands and resources, competing and cooperating on a virtual market. Agents of the basic types are the agents of orders (works), tasks (operations), resources (renewable resources), products (nonrenewable resources), as well as the scene agent. Each agent has its own objective function named the satisfaction function, which is a weighted sum of components that meet various criteria. The value of the system's objective function is refined through the normalized sum of the agents' objective functions. The scheduling algorithm given in [7] includes a negotiation stage used by agents to build a set of conflicting orders for the resources and a conflict resolution stage, which recursively search for placement options taking into account existing limitations. The main advantages of the agent-based scheduling approach proposed in [7] are the agent's knowledge base of preferences when assigning resources to operations and the ability to reschedule in real time if the list of available resources or works has changed.
Heuristic scheduling algorithms based on a genetic algorithm are presented in [8,9]. The genetic algorithm (GA) is a popular optimization technique applicable to various application fields proposed by Goldberg [17]. In [8,9], a sequence of works is encoded into a chromosome, a population of the chromosomes is formed to represent a project, and then the population is transformed via genetic operators. Each of the investigations proposes a fitness function of the total project duration that is to be minimized. In [8], the authors introduce a dense gene concept, which is a fixed chromosome section that encodes a sequence of works that uses the available renewable resources in the most optimal way, providing the smallest remainder of free resources. The scarcity of the renewable resource is determined by solving the simplified scheduling problem with assumption that all the resources are non-renewable and calculating the remaining free resources. In [9], a scheduling problem of prefabricated buildings is considered. The authors consider both renewable resources and prefabricated blocks, or non-renewable resources, and their supplies. The GA consists of two stages: first, a schedule is formed with the assumption that the non-renewable resources are in shortage; second, the non-renewable resource supply's time constraints are added to the schedule and a new search is performed. The advantage of the algorithm presented in [9] is the use of heuristics to generate an initial population that allows a reduction in the number of unfeasible solutions produced by GA.
OptQuest optimizer is a software that incorporates a combination of metaheuristics, such as scatter search, tabu search, and neural networks [18]. As a part of the AnyLogic modeling system [19], OptQuest is used to solve the multi-objective scheduling problem in [10,11]. A simulation model is used as an objective function for the optimization, which in return determines an optimal configuration of input parameters for the simulation model. Two algorithms of the renewable resources' allocation based on dispatching rules or on agent negotiation are compared using the OptQuest optimizer in [10]. Controlled variables are the assignment of resources to works, and the objective function is a multi-object function that estimates the total delay of the works and the total consumption of the renewable resources. The authors concluded that the use of the multi-agent approach decreases the resource waiting time for works but does not increase the throughput. The advantage of the study [11] is that energy consumption is considered as a part of the objective function. The author searches for an effective batch size using OptQuest optimizer and then improves the resources' utilization within the obtained schedule. It is concluded that batch size has a greater impact on the total delay rather than on the energy consumption.
A GA-based solver is embedded into the Tecnomatix Plant Simulation software package [20]. An application of this solver to the scheduling problem is presented in [12,13]. In [12], a flow shop problem is considered with a multi-objective function minimizing the mean flow time and total setup time. The solver optimizes the batch size and sequence of the product batches entering the system. In [13], a real-time job shop scheduling problem is considered with restricted renewable and non-renewable resources and works' due dates. A simulation model is used to assess the idling of processing equipment and the efficiency of the workshop in terms of throughput. The GA solver is used to allocate the renewable and non-renewable resources to the works. Advantages of the study are consideration of both types of resources and the integration of the scheduling phase into a cloud manufacturing system.
Subcontracted renewable resources are considered in [14]. The authors optimize the schedule cost by allocation of own and subcontracted resources taking into account a time-cost contradiction. The method is discussed in detail in Section 5.
We consider the scheduling problem optimizing not only own renewable resources but also subcontracted ones. The problem's objective function is cost minimization of the engaged subcontracted resources. The problem restrictions include the time frame of the earliest and latest work start date and a limited amount of the non-renewable resources available.
The goal of the present paper is to develop and assess an imitation and heuristic (IH) method for the scheduling problem with subcontracted resources. We considered two scheduling theory methods of scheduling on parallel machines by V. S. Tanaev and Y. A. Mezentsev to be most suitable for solving the given problem. We also considered application of a commercial solver to the given problem. Use of simulation and heuristics together allows us to overcome the revealed disadvantages of the considered methods. The developed IH method was applied for schedule search in a project company.

Indices
We consider the following indices: • g is an index of the renewable resource of the allocated competence, g = 1, Q r ; • i is an index of the operation or order contained in the project p, i = 1, N p ; • j is an index of the operation or order that cannot be executed at the same time with the order i because they are to be processed by a single machine g, j = 1, N p ; • k is an index of the time interval specified by the time constraints in the V.C. Tanaev L is an index of the order that can be processed during E k time interval in the V.C. Tanaev method, L = 1, n(k); • p is an index of the project contained in the portfolio, p = 1, P; • r is an index of the renewable resource's competence, r = 1, R; • st is an index of the dynamic programming stage of the Y.A. Mezentsev algorithm, st = 1, N p ; • t is an index of the day, t = τ, T; • v is an index of the non-renewable resource, v = 1, V; • η is an index of the iteration stage of the IH algorithm, η = 1, Ψ.

Sets
We consider the following sets: • e 1 < e 2 < . . . < e k+1 is a set of the variables τ 0 i and dl i values in the V.C. Tanaev method; • N k = i k,1 , i k,2 , . . . , i k,n(k) is a set of orders that can be processed during E k interval with index L in the V.C. Tanaev method; • Pr p is a set of operations in the project p; • Rown is a set of own renewable resources; • Rsc is a set of subcontracted renewable resources; • Z(t, r) is a set of indices for operations performed at time t ≥ τ using own or subcontracted renewable resource r; • Z(t, v) is a set of indices for operations performed at time t ≥ τ using non-renewable resource v; • Θ(t, r) is a set of indices for operations performed at time t ≥ τ using own resource r; • Y(t, r) is a set of indices for operations performed at time t ≥ τ using subcontracted renewable resource r.

Parameters
We identified the following parameters: • c p,i,r is a duration of the time interval M p,i,r (η) in the IH algorithm; • d i > 0 is an operation duration in case of one project in the portfolio, i.e., P = 1; • dl i > 0 is a deadline for processing the order i; it can be calculated as: Q r is the available quantity of the own renewable resource with competence r; • Q t,v is the current quantity of each non-renewable resource v at time t; • q p,i,r ≥ 0 is the required amount of renewable resources of allocated competence r to perform the operation i of the project p; q p,i,r = 0 when the renewable resource r is not required to process the operation i of the project p; • q − p,i,v ≤ 0 and q + p,i,v ≥ 0 are the amounts of non-renewable resource v consumed when the operation i of the project p starts and produced when the operation i ends, respectively; q − p,i,v = 0 when the non-renewable resource v is not required to start the operation i of the project p; the variable q + p,i,v = 0 when the non-renewable resource v is not produced as an outcome of the operation i of the project p; • SC p,t,r is a cost of the subcontracted resource of allocated competence r engaged in the project p at the time t; • sr i,r is the daily cost of performing the operation i using a unit quantity of the renewable resource r in case of one project in the portfolio, i.e., P = 1; • sr p,i,r ≥ 0 is the daily cost of performing the operation i of the project p using a unit quantity of the renewable resource of allocated competence r; • sr p,i,r ≥ 0 is the total cost of performing the operation i of the project p using a subcontracted renewable resource r; • Tdl is the global project's deadline that must not be exceeded, Tdl = max i dl i ; • U p,t,r is the percentage utilization of the renewable resource r engaged in the project p at the time t.
is the average utilization of the renewable resource r for the operation i of the project p during the time interval δ p,i,r (η) in the IH algorithm; is the average utilization of the renewable resource r for the operation i of the project p during the time interval δ − p,i,r (η) in the IH algorithm; is the average utilization of the renewable resource r for the operation i of the project p during the time interval δ + p,i,r (η) in the IH algorithm; shifted to the left on the time axis in the IH algorithm; shifted to the right on the time axis in the IH algorithm; • λ g,i is the cost of processing the order i by the resource g; is a function that shows the presence of subcontracted resources g at the time t, g = Q r , M; • ξ p,i,r is a threshold of the sr p,i,r total cost of performing the operation i of the project p using a subcontracted renewable resource r. • σ(t) is a function of orders allocated to machines or renewable resources in the V.C. Tanaev method; • τ 0 p,i ≥ τ and τ 1 p,i ≥ τ are the earliest and latest possible start times given for the operation i of the project p; τ 0 p,i = τ 1 p,i = τ if the operation i is only allowed to start at the time τ; • τ 0 i and τ 1 i are orders' earliest and latest possible start times in case of one project in the portfolio, P = 1; •τ g,i is the actual delay between start of the i-th order on the g-th machine upon the previous order completion in the Y.A. Mezentsev algorithm; st is the completion time by the resource g for all the orders that exist at the stages from the first to the st-th in the Y.A. Mezentsev algorithm; • ω is the cost of processing the order per time unit in conventional units; • Ψ is the predefined maximum number of the IH algorithm steps.

Decision Variables
The decision variables may vary regarding the different scheduling problem statements. We identify the following decision variables connected with the problem: are integer variables of a sequence of orders to perform by each renewable resource g at each time interval E k for the scheduling problem given in Section 5. executed by the machine g; otherwise, it equals 0. In Section 5.3, we assume that the machine g can be an own or subcontracted one; • ε i,j are binary variables for the scheduling problem given in Section 5.3. The variable indicates orders that cannot be executed at the same time because they are to be processed by a single machine g. The variable ε i,j equals 1 if the order i is to be completed before the order j; otherwise, it equals 0; (i, j) ∈ Pr p .

Scheduling Problem Statement
The scheduling problem statement is given by the authors in [21]. Below are the key points. The notations are presented in Section 3.
We assume that a set of operations in each project p appears in order of increasing of the operation's cost sr p,i,r ; sr p,i,r = 0 if the subcontracted resource r is not required by the operation i of the project p.
We denote the set of indices for the operations performed at the time t ≥ τ utilizing the renewable resource r as follows: We denote the set of indices for the operations performed at the time t ≥ τ utilizing the non-renewable resource v as follows: We denote the set of indices for the operations performed at the time t ≥ τ utilizing the own renewable resource r within the available amount Q r : The set of indices for the operations performed at the time t ≥ τ utilizing the subcontract resource r is defined as follows: The percentage utilization of the resource r engaged in the project p at the time t is defined as follows: The cost of the subcontract resource of an allocated competence r involved in the project p at the time t is defined by the formula: The current volume of the non-renewable resource v at the time t is defined as: The scheduling problem can be formalized as follows: The objective function (8) minimizes the total cost of the subcontracted resources in the case of exceeding the availability of the own ones. The constraint (9) ensures availability of the required amount of non-renewable resources at the time of the operation start. The constraint (10) imposes a time frame on the start dates of the operations.

Scheduling Methods Application
We consider the schedule optimization problem for the parallel system having identical machines in the presence of delays processing the orders. The problem is one of the closest problems studied by the scheduling theory. According to [15], the problem is formalized as follows: Pm/r j /C max . For the problem, a parametric algorithm of dynamic programming with an alternative dropout option was proposed by Y. A. Mezentsev et al. [4].
Another closest problem studied by the scheduling theory is a parallel system schedule optimization using identical machines in the presence of orders' due dates. The problem is formalized as follows: Pm/brkdwn, r j / ∑ ω j U j . As a solution to the problem, a schedule construction algorithm using the given due dates is proposed by V. S. Tanaev et al. [5].
We also considered Lingo commercial solver's application to the scheduling problem according to the scheme given in [14]. Here, the problem is modelled as a mixed binary linear program that minimizes the project cost, including subcontracted and own resources cost.
We applied the considered algorithms to solve a trial small-scale scheduling problem using two own renewable resources Q 1 = 2 of the same competence R = 1, and one project P = 1 with the number of the operations N 1 = 7. Seven orders, or operations, are fed at a different time into a parallel system with two identical machines, or renewable resources. Table 1 contains the initial information about the operation's duration and time frame between order processing's earliest possible start and latest finish. We assume the operations are ordered by the earliest start time. The progress of solving the problem is given below.

Method Developed by Y. A. Mezentsev
We applied a parametric dynamic programming algorithm with the optional exclusion of the found alternatives to the considered scheduling problem. The notations are given in Section 3. The decision variables are Boolean variables y g,i allocating the order i to the machine g.
The dynamic algorithm variables can be calculated as given in the study [4]: Mathematics 2021, 9, 2098 9 of 22 The parametric algorithm includes the following stages: 1.
Input of the initial data τ 0 i , d i , dl i , i = 1, N p and H, K , Q r , and N p parameters. Set If st > N p then go to Point 7; 4.
Generate all the feasible schedules and calculate f g,st τ 0 st , d st , y g,st and the schedule length Check the number of the generated schedules N st . If N st ≤ K then go to Point 2; otherwise, go to Point 6; 6.
Discard Q r H−1 out of the schedules generated at Point 4 with the maximum schedule length ϕ g,st τ 0 i , d i , y g,i . Go to Point 2; 7.
Choose the schedules with the minimum makespan. Determine the calendar plan by reverse dynamic programming.
We consider application of the parametric algorithm to the trial scheduling problem given in Table 1. Examples of the calculated algorithm characteristics are given in Tables 2  and 3 for the first and last algorithm stages. We set H = 2, K = 2 2 = 4. Table 2. Results of the first algorithm stage completion. Table 3 contains the dark filled cells connected to the schedules providing the minimum local makespan on the given stage. Table 4 contains the four schedules identified as the solutions of the given small-scale scheduling problem.
Since the machines are identical, schedule 1 is equal to schedule 4 and schedule 2 is equal to schedule 3. Schedules 1 and 2 differ by two last orders allocated to the opposite machines. All of the schedules reveal deadline violation on orders 6 and 7.
The obtained schedules are shown in Figure 1 in the form of a Gantt chart. As we can see, the schedules found by the Y. A. Mezentsev method for a parallel system with identical machines and orders' earliest start time do not support the orders' due date. Subcontracted resources are not considered by the method. Thus, = 0 for the schedule found by Y. A. Mezentsev.

Method Developed by V. C. Tanaev
We considered an algorithm implemented by V. C. Tanaev. The notations used are given in Section 3.
The decision variables are schedules , ( ) for each resource at each time interval . These are integer type variables.
Two cases can occur: 1) ≤ if the existing machines are enough to process all the orders within time intervals ; 2) > if the existing machines are not enough to process all the orders at the time intervals and it is necessary to attract an amount ≤ ( − ) of subcontracted machines. An algorithm is given in [5] and contains the following stages. We modified the algorithm by adding the stages 4 to 7 to reduce the number of interruptions:

1.
= 0; As we can see, the schedules found by the Y. A. Mezentsev method for a parallel system with identical machines and orders' earliest start time do not support the orders' due date. Subcontracted resources are not considered by the method. Thus, Ft = 0 for the schedule found by Y. A. Mezentsev.

Method Developed by V. C. Tanaev
We considered an algorithm implemented by V. C. Tanaev. The notations used are given in Section 3.
The decision variables are schedules S k,g (t) for each resource g at each time interval E k . These are integer type variables.
Two cases can occur: (1) M ≤ Q r if the existing machines are enough to process all the orders within time intervals E k ; (2) M > Q r if the existing machines are not enough to process all the orders at the E k time intervals and it is necessary to attract an amount M ≤ (M − Q r ) of subcontracted machines.
An algorithm is given in [5] and contains the following stages. We modified the algorithm by adding the stages 4 to 7 to reduce the number of interruptions: The next index of the time interval: k = k + 1. If k > β then go to Point 9; otherwise, perform the following actions. Form a set of orders N k = i k,1 , i k,2 , . . . , i k,n(k) , which can be processed at the E k interval. The orders are sorted by descending the cost of the processing per time unit; 3.
If min i∈N k {d i } < ∆ k then go to Point 8; otherwise, go to Point 4;

4.
Start to browse a set of machines with given competence: g = 0; 5.
The next index of the machine: g = g + 1. If (g > M) OR (N k = ∅), then go to Point 2; otherwise, start to browse the set of orders from the beginning N k : L = 0; 6. The next index of the order: L = L + 1. If L > n(k), then go to Point 5; otherwise, go to Point 7; 7.
Assign the machine g to the order i k,L ∈ N k and form the schedule S k,g (t) at the E k time interval: a. If (k = 1) OR (k = 1 AND the order i k,L at the previous time interval E k−1 has not been assigned to the machine g = g) OR (k = 1 AND the order i k,L at the previous time interval E k−1 has been assigned to the machine g = g, g > Q r AND n(k) ≤ Q r ), then assign the machine g to the order i k,L : i.
Eliminate the order i k,L from the set N k . Go to Point 5; b. If (k = 1 AND the order i k,L at the previous time interval E k−1 has been assigned to the machine g = g, g ≤ Q r ) OR (k = 1 AND the order i k,L at the previous time interval E k−1 has been assigned to the machine g = g and g > Q r AND n(k) > Q r ), then go to Point 6;

8.
Calculate the processing duration of all the orders constituting the N k set according to the formula: Apply a packing algorithm intended to assign the N k set orders having d i durations to the set M machines for the E k time interval. At this stage, we consider the following statements are true: a. Calculate at the time interval (e k ; e k + M∆ k ] a function σ(t), where: Form the schedule S k,g (t) = σ(t + (g − 1)∆ k ), g = 1, M; c.
The end of the algorithm.
We used the algorithm to solve the scheduling problem given in Table 1. We sorted the set N p orders according to the information given in Table 1.
2. The next index of the time interval: = + 1. If > then go to Point 9; otherwise, perform the following actions. Form a set of orders = , , , , … , , ( ) , which can be processed at the interval. The orders are sorted by descending the cost of the processing per time unit; 3. If min ∈ { } < ∆ then go to Point 8; otherwise, go to Point 4; 4. Start to browse a set of machines with given competence: = 0; 5. The next index of the machine: = + 1. If ( > ) OR ( = ∅), then go to Point 2; otherwise, start to browse the set of orders from the beginning : = 0; 6. The next index of the order: = + 1. If > ( ), then go to Point 5; otherwise, go to Point 7; 7. Assign the machine to the order , ∈ and form the schedule , ( ) at the time interval: a. If ( = 1) OR ( ≠ 1 AND the order , at the previous time interval has not been assigned to the machine ≠ ) OR ( ≠ 1 AND the order , at the previous time interval has been assigned to the machine ≠ , > AND ( ) ≤ ), then assign the machine to the order , : i.
= − ∆ , , ( ) = , , where = ( ; ]; ii. Eliminate the order , from the set . Go to Point 5; b. If ( ≠ 1 AND the order , at the previous time interval has been assigned to the machine ≠ , ≤ ) OR ( ≠ 1 AND the order , at the previous time interval has been assigned to the machine 9. The end of the algorithm.
We used the algorithm to solve the scheduling problem given in Table 1. We sorted the set orders according to the information given in Table 1.  Results of the method's execution upon the trial problem are given in Table 5.
Since max ( ) = 4, the number of machines is defined as four. However, the maximum number of the machines used simultaneously turned out to be three as two orders were assigned to the machine at the time interval during schedule searching. Hence, we set = 3, = 2, and = 1. Results of the method's execution upon the trial problem are given in Table 5.
Since max k n(k) = 4, the number of machines is defined as four. However, the maximum number of the machines used simultaneously turned out to be three as two orders were assigned to the machine r 1 at the time interval E 6 during schedule searching. Hence, we set M = 3, Q r = 2, and M = 1. The schedule formed by the V. S. Tanaev method is shown in Figure 3.

Interval and Its
For the schedule found by the V. S. Tanaev method, the cost of the subcontracted resources is Ft = 4ω.
The schedule found is not optimal according to the objective function (8). An example of a more effective schedule satisfying all the requirements (9)-(10) is shown in Figure 4; the subcontracted resources cost is Ft = 2ω.
For the schedule found by the V. S. Tanaev method, the cost of the subcontracted resources is = 4 .
The schedule found is not optimal according to the objective function (8). An example of a more effective schedule satisfying all the requirements (9)-(10) is shown in Figure 4; the subcontracted resources cost is = 2 . Without loss of generality, let us identify = 100. Then, = 400 for the schedule found by the V. C. Tanaev method.

Method Based on the Commercial Solver Application
We applied a method based on a Lingo solver [22] to the trial scheduling problem. We reduced our mathematical model into a model suitable for linear programming [14].
The decision variables are the following: • are integer variables indicating the start times of the orders, ∈ ; • , are binary variables indicating what renewable resource , either own or subcontracted, is allocated to a particular order ; • , are binary variables indicating orders that cannot be executed at the same time because they are to be processed by a single resource .

Method Based on the Commercial Solver Application
We applied a method based on a Lingo solver [22] to the trial scheduling problem. We reduced our mathematical model into a model suitable for linear programming [14].
The decision variables are the following: • x i are integer variables indicating the start times of the orders, i ∈ Pr p ; • y g,i are binary variables indicating what renewable resource g, either own or subcontracted, is allocated to a particular order i; • ε i,j are binary variables indicating orders that cannot be executed at the same time because they are to be processed by a single resource g.
The mathematical model of this problem applied to a project is described as follows: ∑ g∈Rown∪R sc x i + d i ≤ x j + ML · 1 − ε i,j + ML · 2 − y g,i − y g,j , ∀(i, j) ∈ Pr p , ∀g ∈ Rown ∪ Rsc, (17) For the problem given in Table 1, the following variable values are set: • The number of projects is P = 1; • The cardinality of the Rown set is |Rown| = 2; • The cardinality of the Rsc set is |Rsc| = 7, where 7 is the number of orders that can be processed simultaneously; • The cardinality of the Pr p set is Pr p = 7; • Ksc = 400 that is limited by the maximum value of the subcontracted resources cost found by the V. S. Tanaev and Y. A. Mezentzev methods; • ML = 10000; • Tdl = 10 days according to the dl i values from Table 1; • λ g,i = 1, i f g ∈ Rown 100, i f g ∈ Rsc ∀ i ∈ Pr p ; we assume that the operation cost should be much lower when processed by the own resource comparing to the subcontracted.
It is defined for this problem that limited non-renewable resources are not considered while they are being introduced in the Formulation (8)- (10).
As a result of Lingo execution upon the linear programming problem (15)-(23), a set of alternative optimal schedules was formed providing the same minimum total used resources cost Kres = 7. Two optimal schedules are presented in Figure 5.
For the problem given in Table 1, the following variable values are set: • The number of projects is = 1; • The cardinality of the set is | | = 2; • The cardinality of the set is | | = 7, where 7 is the number of orders that can be processed simultaneously; • The cardinality of the set is = 7; • = 400 that is limited by the maximum value of the subcontracted resources cost found by the V. S. Tanaev and Y. A. Mezentzev methods; • = 10000; • = 10 days according to the values from Table 1; ; we assume that the operation cost should be much lower when processed by the own resource comparing to the subcontracted.
It is defined for this problem that limited non-renewable resources are not considered while they are being introduced in the Formulation (8)- (10).
As a result of Lingo execution upon the linear programming problem (15)-(23), a set of alternative optimal schedules was formed providing the same minimum total used resources cost = 7. Two optimal schedules are presented in Figure 5. As we can see, the commercial solver was able to reach subcontracting cost = 0 for the schedules found, but there are violations of the deadline and the orders' earliest start time. Thus, application of the commercial solver method leads to violation of the restrictions (9) and (10) of the problem considered.

Imitation and Heuristic Method
A key concept of the imitation and heuristic algorithm is integration of processes' imitation and some heuristic rules to improve the initial schedule. The IH method algorithm is based on an application of a multi-agent resource conversion process (MRCP) model [23]. The MRCP model is intended to describe discrete processes converting input non-renewable resources into output ones using renewable resources, or machines, throughout a given time interval.
An agent of the MRCP model is a decision maker model having formalized knowledge about resources' allocation using production rules. The MRCP model also includes a logistics agent. The logistics agent controls the current value and lifetime of the non-renewable resources and ensures fulfillment of the restriction (9) by launching the purchase or production process of the non-renewable resource required in case its current volume is decreased to a critical value or the resource's lifetime is exceeded. Figure 5. Two optimal schedules provided by the commercial solver.
As we can see, the commercial solver was able to reach subcontracting cost Ft = 0 for the schedules found, but there are violations of the deadline and the orders' earliest start time.
Thus, application of the commercial solver method leads to violation of the restrictions (9) and (10) of the problem considered.

Imitation and Heuristic Method
A key concept of the imitation and heuristic algorithm is integration of processes' imitation and some heuristic rules to improve the initial schedule. The IH method algorithm is based on an application of a multi-agent resource conversion process (MRCP) model [23]. The MRCP model is intended to describe discrete processes converting input non-renewable resources into output ones using renewable resources, or machines, throughout a given time interval.
An agent of the MRCP model is a decision maker model having formalized knowledge about resources' allocation using production rules. The MRCP model also includes a logistics agent. The logistics agent controls the current value and lifetime of the nonrenewable resources and ensures fulfillment of the restriction (9) by launching the purchase or production process of the non-renewable resource required in case its current volume is decreased to a critical value or the resource's lifetime is exceeded.
The IH algorithm is a cycle with alternating stages of imitation and application of heuristic rules to improve the schedule. During the imitation stage, the schedule is fed to the MRCP model input, and the model evaluates the subcontracted resources cost according to the Formula (8). During the heuristic stage, the algorithm shifts the start days for the operations, where the subcontracted resources cost exceeds the given threshold. During shifting, the restriction (10) is ensured. The algorithm stops either in the absence of exceeding the threshold of the subcontracted resources cost, or in case a certain number of cycles is reached.
We assess the number of iterations of the IH algorithm. The number Ψ of search steps is a parameter of the algorithm. For each step, the algorithm sequentially bypasses all operations of the project portfolio to identify bottlenecks and eliminate them.
The number of operations of the project portfolio is calculated as follows: N = ∑ P p=1 N p . Thus, the complexity of the proposed IH algorithm linearly depends on the number of operations in the project portfolio according to the formula: Ψ·N. The quality of the schedule found is defined based on the value of the function (8) when restrictions (9)-(10) are fulfilled.

Case Study Results
The modified IH scheduling method was implemented in the BPsim software package including the BPsim simulation system and BPsim decision support system [24] with a common database. The BPsim simulation system is used to develop and simulate the MRCP model of the processes under investigation. The system supports graphical MRCP notation, where a user can identify the types of nodes, i.e., operations and agents, logical links between them, and list the available non-renewable and renewable resources.
The following parameters can be assigned to each operation: the duration, start condition, amount of non-renewable and renewable resources required to complete the operation, and amount of produced non-renewable resources. To each agent, the behavior rules can be described in a form of if-then rules and assigned. The agents can affect the amount of resources and operations' start conditions. The BPsim decision support system includes decision search diagrams based on the same resources used to build the simulation model. Development of the diagrams is based on the UML sequence diagrams [25] and Transact-SQL database management language [26]. The decision search diagrams are used to provide a visual comparison for multiple alternative decisions utilizing implemented user rules. The diagrams were applied to implement the algorithm of the IH scheduling method.
A case study was used to assess the subcontracted resources cost depending on the current schedule of 10 projects with 35 operations in a project company. The company has its own renewable resources, i.e., a staff of eight people with different skills (competencies). The following competencies of the staff were defined: documentation design (three people), carrying out an installation work (four people), and material supply (one person). The nonrenewable resources of the company are the construction objects. The operation duration varied from 6 days to 90 days depending on the type of operation. The scheduling time interval was 430 days.
Following the proposed IH scheduling method, a simulation model was developed in the BPsim simulation system. The model's inputs are the labor costs and the agreed start dates for each operation, as well as the schedule. Assessment of the subcontracted resources cost was produced by the model.
For the initial schedule, the model outcome is presented in Figure 6a. For the schedule provided by the IH method, the model outcome is presented in Figure 6b.
The figures contain the percentage utilization of the three own resources with different competence marked with the red, blue, and green lines. When the blue and red resources utilization reaches the 100% level, it means the subcontracted resources with the matched competence are used during this time. The green resource utilization does not exceed 30% in both figures, so subcontracted resources with the same competence are not required.
The figures contain the percentage utilization of the three own resources with different competence marked with the red, blue, and green lines. When the blue and red resources utilization reaches the 100% level, it means the subcontracted resources with the matched competence are used during this time. The green resource utilization does not exceed 30% in both figures, so subcontracted resources with the same competence are not required.
(a) (b) Figure 6. Percentage utilization U of own resources iRes performing the project's portfolio depending on time: (a) Initial schedule; (b) Schedule provided by the IH method. Here, iRes2 is a resource with the documentation design competence, iRes3 is a resource with the competence to carry out installation works, and iRes4 is a resource with material supply competence.
The initial and the final schedules coincide up to 110 days due to the low resource utilization at the initial schedule during these days. Therefore, the IH algorithm has not shifted the operations' start day during the first 110 days in the final schedule. Table 6 contains allocation of the subcontracted resources' cost in rubles and man per day (m/d) on the project's operations for the initial schedule. For the initial schedule, the subcontracted resources are required for the projects 1, 3, 4, 5, and 9. The IH method allowed us to highlight several operations with a subcontracted resources cost exceeding the critical value equal to 60,000 rubles. The ones are installation works for projects 3 and 9 (95,990 and 70,953 rubles accordingly), and all operations for project 4 (65,998 rubles). Table 7 contains allocation of the subcontracted resources' cost in rubles and m/d on the project's operations for the schedule provided by the IH method.
The IH method is proposed resource allocation and shifting of the operations' start date within a given timeframe. Although, the makespan of the found schedule was increased compared with the initial schedule, all operations and project deadlines were met. Here, iRes2 is a resource with the documentation design competence, iRes3 is a resource with the competence to carry out installation works, and iRes4 is a resource with material supply competence.
The initial and the final schedules coincide up to 110 days due to the low resource utilization at the initial schedule during these days. Therefore, the IH algorithm has not shifted the operations' start day during the first 110 days in the final schedule. Table 6 contains allocation of the subcontracted resources' cost in rubles and man per day (m/d) on the project's operations for the initial schedule. For the initial schedule, the subcontracted resources are required for the projects 1, 3, 4, 5, and 9. The IH method allowed us to highlight several operations with a subcontracted resources cost exceeding the critical value equal to 60,000 rubles. The ones are installation works for projects 3 and 9 (95,990 and 70,953 rubles accordingly), and all operations for project 4 (65,998 rubles). Table 7 contains allocation of the subcontracted resources' cost in rubles and m/d on the project's operations for the schedule provided by the IH method.
The IH method is proposed resource allocation and shifting of the operations' start date within a given timeframe. Although, the makespan of the found schedule was increased compared with the initial schedule, all operations and project deadlines were met.
It should be noted that the objective function considered is related to minimization of the subcontracted resources' cost while meeting restrictions. As a result, the subcontracted resources' cost of the schedule found was reduced compared with the initial schedule in terms of man per day from 42 m/d to 9 m/d and in terms of rubles from 249,936 rubles to 53,453 rubles for half a year, i.e., more than four times for both outcomes. For the scheduling problem with 35 operations, eight renewable resources and one type of non-renewable resource, the running time of the HI algorithm is estimated at 6 min. This time is comparable to the running time of the methods based on agents, OptQuest, and Plant Simulation application. The HI algorithm running time is much more than the one of the methods proposed by Y. A. Mezentsev and V. S. Tanaev as well as the Lingo-and GA-based ones, estimated at tens of seconds. It should be noted that the search time for the simulation-based methods is always more than the search time for the methods based on the optimization algorithms. When simulating, a time for conducting one experiment varies from a few seconds to several minutes depending on the model dimension and simulation tool used. In case of solving the optimization problem, it is necessary to search for solutions in a search space while the alternative solution is estimated by conducting an experiment with the model; therefore, the total search time increases to tens or even hundreds of minutes. Nevertheless, simulation-based methods have the advantage of being able to determine the objective function and the constraints required without a reduction to a specific mathematical model used in optimization algorithms.
We also applied the HI algorithm to the scheduling problem with a large data set. A construction holding was considered with 302 building operations, 119 renewable resources or construction machinery, and 161 types of non-renewable resources or construction material. The detailed holding description is given in [27]. For the problem with a large data set, the running time of the HI algorithm is estimated at one hour and 38 min. The given duration of the algorithm operation is rather long compared to the obtained one for the small dimension problem but is acceptable for decision making that is not in real time.

Conclusions
The scheduling problem was formulated with the objective function of minimization of the subcontracted resources cost, presence of restrictions on non-renewable resources, and operations' earliest start time and due dates. Analysis of the different scheduling methods based on the scheduling theory, optimization algorithms, and agent-based simulation was conducted. The analysis revealed factors preventing application of the methods to the problem under consideration. The factors include lack of a search for the optimal allocation of renewable resources on operations in terms of minimizing the subcontracted resources' cost and lack of accounting for limited non-renewable resources. The given factors indicate the relevance of the development of the hybrid scheduling method based on integration of simulation and heuristic modeling.
The IH method algorithm was developed based on the multi-agent simulation model of the resource's allocation with operations performing and heuristic rules of the operations start days shifting. The IH method is used to search the schedule that meets the time and resources restrictions and has a minimum subcontracted resources cost.
A case study was conducted to assess the subcontracted resources cost for the real scheduling problem of a project company. Application of the IH method allowed us to compose a schedule that reduced the company's waste on subcontracted resources by more than four-fold for half a year compared with the schedule provided by decision makers based on their knowledge and experience.
Comparison of the IH method and the other scheduling heuristic methods was performed. The conditions were identified, under which the new IH method is more effective than the other ones. The conditions include a focus on optimizing the project portfolio cost with a fixed portfolio duration. In future, the authors plan to refine the IH method for solving the multicriteria problem of finding a schedule that is optimal in terms of the makespan and the renewable resources cost.