*Article* **A Static Assignment Algorithm of Uniform Jobs to Workers in a User-PC Computing System Using Simultaneous Linear Equations**

**Xudong Zhou 1, Nobuo Funabiki 1,\*, Hein Htet 1, Ariel Kamoyedji 1, Irin Tri Anggraini 1, Yuanzhi Huo <sup>1</sup> and Yan Watequlis Syaifudin <sup>2</sup>**


**\*** Correspondence: funabiki@okayama-u.ac.jp

**Abstract:** Currently, the *User-PC computingsystem (UPC)* has been studied as a low-cost and highperformance distributed computing platform. It uses idling resources of personal computers (PCs) in a group. The job-worker assignment for minimizing *makespan* is critical to determine the performance of the UPC system. Some applications need to execute a lot of *uniform jobs* that use the identical program but with slightly different data, where they take the similar CPU time on a PC. Then, the total CPU time of a worker is almost linear to the number of assigned jobs. In this paper, we propose a *static assignment algorithm* of *uniform jobs* to workers in the UPC system, using *simultaneous linear equations* to find the *lower bound* on *makespan*, where every worker requires the same CPU time to complete the assigned jobs. For the evaluations of the proposal, we consider the uniform jobs in three applications. In *OpenPose*, the CNN-based *keypoint* estimation program runs with various images of human bodies. In *OpenFOAM*, the physics simulation program runs with various parameter sets. In *code testing*, two open-source programs run with various source codes from students for the *Android programming learning assistance system (APLAS)*. Using the proposal, we assigned the jobs to six workers in the testbed UPC system and measured the CPU time. The results show that *makespan* was reduced by 10% on average, which confirms the effectiveness of the proposal.

**Keywords:** UPC; distributed computing platform; uniform job; static assignment; linear equations

#### **1. Introduction**

Currently, the *User-PC computing (UPC)* system has been studied as a low-cost and high-performance distributed computing platform [1]. The UPC system uses idling resources of personal computers (PCs) in a group to handle a number of various computing jobs from users. Then, the proper assignment of incoming jobs to workers is very important to effectively deal with them by using computational resources properly. As a result, the job assignment algorithm is critical to achieve the minimization for *makespan* to complete all the demanded jobs in the UPC system.

Previously, we proposed the algorithm of assigning *non-uniform jobs* to workers in the UPC system [2]. In *non-uniform jobs*, the programs are much different from each other, including the developed programming languages, the number of threads, and the requiring data. The execution time for each *non-uniform job* is highly different from the others. The previous algorithm can find the job-worker assignment through two stages sequentially, of which are heuristic due to the nature of the *NP-hardness* and cannot guarantee the optimality of the solution.

Some applications need to execute a lot of *uniform jobs* that use the identical program but with slightly different data/files, where they take a similar CPU time on a PC. The applications include deep learning (machine learning), physics simulations, software testing, computer network simulations, mathematical modeling, and mechanics modeling.

**Citation:** Zhou, X.; Funabiki, N.; Htet, H.; Kamoyedji, A.; Anggraini, I.T.; Huo, Y.; Syaifudin, Y.W. A Static Assignment Algorithm of Uniform Jobs to Workers in a User-PC Computing System Using Simultaneous Linear Equations. *Algorithms* **2022**, *15*, 369. https:// doi.org/10.3390/a15100369

Academic Editor: Frank Werner

Received: 28 August 2022 Accepted: 4 October 2022 Published: 7 October 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

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

These jobs have the common feature of a similar CPU time when they run on a specific PC. The *uniform jobs* often need a long CPU time. For example, in physics or network simulations, it can take several days to run one job. Nevertheless, it will be necessary to find the best result of all the input data by repeating them to slightly change some parameter values for the program and running them. This work can be common in research activities using computer simulations.

In this paper, we propose a *static assignment algorithm* of *uniform jobs* to workers in the UPC system, using *simultaneous linear equations* to find the *lower bound* on *makespan*, where every worker requires the same CPU time to complete the assigned jobs. The *simultaneous linear equations* describe the equality of the estimated CPU time among the workers, and the equality of the total number of assigned jobs to workers with the number of given jobs. The estimated CPU time considers simultaneous executions of multiple jobs on one worker by using its multiple cores. Since solutions of *simultaneous linear equations* become real numbers in general, the integer number of jobs assigned to each worker is introduced to them in a greedy way.

For evaluations of the proposal, we consider *uniform jobs* in the three applications for the UPC system, namely, *OpenPose* [3], *OpenFOAM* [4], and *code testing* [5,6]. For *OpenPose*, the CNN-based program runs with 41 images of human bodies. For *OpenFOAM*, the physics simulation program runs with 32 parameter sets. For *unit testing*, the opensource programs run with 578 source codes that were submitted from students to the server in the *Android programming learning assistance system (APLAS)*. These jobs were applied to the proposed algorithm and were assigned to six workers in the testbed UPC system by following the results. Then, the CPU time was measured by running them. For comparisons, two simple algorithms were also implemented where the jobs were applied, and the CPU time was measured. The evaluation results show that the difference between the longest CPU time and the shortest one among the six workers became 92 s, and *makespan* of the UPC system was reduced by 10% on average from the results by comparative algorithms. Thus, the effectiveness of the proposal was confirmed.

The proposed algorithm limits the application to the jobs where the CPU time is nearly equal to a worker. This limitation can simplify the job scheduling algorithm to only considering the number of jobs assigned to each worker, while neglecting the differences between individual jobs. Fortunately, it is possible to alleviate this limitation to a certain degree by considering the granularity of the CPU time on a worker. The CPU time of a job that is applicable to the proposal is often proportional to the number of iteration steps before the termination, or to the number of elements in the computational model. For example, in computer network simulations, the number of iteration steps need to be selected with the unit time before simulations, where the CPU time is usually proportional to it. By considering a multiple of a constant number of iteration steps, such as 100, the CPU time can be estimated even if the number of iteration steps is widely changed with this granularity. In future works, we will study this extension of the proposed algorithm to increase its applicable applications.

The rest of this paper is organized as follows: Section 2 discusses related works. Section 3 reviews the UPC system, *OpenPose*, *OpenFOAM*, and *code testing* in *APLAS*. Section 4 presents the *static assignment algorithm* of *uniform jobs* to workers in the UPC system. Section 5 evaluates the proposal through experiments. Section 6 extends the proposal to multiple job-type assignments. Finally, Section 7 concludes this paper with future works.

#### **2. Related Works in the Literature**

In this section, we discuss some related works in the literature.

In [7], Lin proposed several linear programming models and algorithms for identical jobs (*uniform jobs*) on parallel uniform machines for individual minimizations of several different performance measures. The proposed linear programming models provide structured insights of the studied problems and provide an easy way to tackle the scheduling problems.

In [8], Mallek et al. addressed the problem of scheduling identical jobs (*uniform jobs*) on a set of parallel uniform machines. The jobs are subjected to conflicting constraints modeled by an undirected graph *G*, in which adjacent jobs are not allowed to be processed on the same machine. The minimization of the maximum *makespan* in the schedule is known to be *NP-hard*. To solve the general case of this problem, they proposed mixed-integer linear programming formulations alongside lower bounds and heuristic approaches.

In [9], Bansal et al. proposed the two-stage *Efficient Refinery Scheduling Algorithm (ERSA)* for distributed computing systems. In the first stage, it assigns a task according to the min–max heuristic. In the second stage, it improves the scheduling by using the refinery scheduling heuristic that balances the loads across the machines and reduces *makespan*.

In [10], Murugesan et al. proposed a multi-source task scheduler to map the tasks to the distributed resources in a cloud. The scheduler has three phases: the task aggregation, the task selection, and the task sequencing. By using the ILP formulation, this scheduler minimizes *makespan* while satisfying the budget allotted by the cloud user based on the divisible load theory.

In [11], Garg et al. proposed the *adaptive workflow scheduling (AWS)* for grid computing using the dynamic resources based on the rescheduling method. The AWS has three stages of the initial static scheduling, the resource monitoring, and the rescheduling, to minimize *makespan* using the directed acyclic graph workflow model for grid computing. It deals with the heterogeneous dynamic grid environment, where the availability of computing nodes and link bandwidths are inevitable due to existences of loads.

In [12], Gawali et al. proposed the two-stage *Standard Deviation-Based Modified Cuckoo Optimization Algorithm (SDMCOA)* for the scheduling of distributed computing systems. In the first stage, it calculates the sample initial population among all the available number of task populations. In the second stage, the modified COA immigrates and lays the tasks.

In [13], Bittencourt et al. reviewed existing scheduling problems in cloud computing and distributed systems. The emergence of distributed systems brought new challenges on scheduling in computer systems, including clusters, grids, and clouds. They defined a taxonomy for task scheduling in cloud computing, namely, pre-cloud schedulers and cloud schedulers, and classified existing scheduling algorithms in the taxonomy. They introduced future directions for scheduling research in cloud computing.

In [14], Attiya et al. presented a modified *Harris hawks optimization (HHO)* algorithm based on the *simulated annealing (SA)* for scheduling the jobs in a cloud environment. In this approach, SA is employed as a local search algorithm to improve the convergence rate and the solution quality generated by the standard HHO algorithm. HHO is a novel populationbased, nature-inspired optimization paradigm proposed by Heidari et al. [15]. The main inspiration of HHO is the cooperative behavior and the chasing style of Harris' hawks in nature. In the HHO model, several hawks explore prey, respectively, and simultaneously after attacking the target from different directions to surprise it.

In [16], Al-Maytami et al. presented a novel scheduling algorithm using *Directed Acyclic Graph (DAG)* based on the *Prediction of Tasks Computation Time algorithm (PTCT)* to estimate the preeminent scheduling algorithm for prominent cloud data. The proposed algorithm provides a significant improvement with respect to *makespan* and reduces the computational complexity via employing *Principal Components Analysis (PCA)* and reducing the *Expected-Time-to-Compute (ETC)* matrix.

In [17], Panda et al. proposed an *energy-efficient task scheduling algorithm (ETSA)* to address the demerits associated with the task consolidation and scheduling. The proposed algorithm *ETSA* takes into account the completion time and the total utilization of a task on the resources, and follows a normalization procedure to make a scheduling decision. The *ETSA* provides an elegant trade-off between energy efficiency and *makespan*, more so than the existing algorithms.

#### **3. Reviews of UPC System and Three Applications**

In this section, we review the *User-PC computing system (UPC)* system and the three applications in this paper.

#### *3.1. UPC System*

First, we review the *UPC* system. The UPC system can provide computational powers efficiently for members in a group, such as engineers in a company or students in a laboratory, by using idling computing resources of their PCs. To allow various application programs to run on different PC environments, the UPC system adopts *Docker*. *Docker* is a popular software tool that has been designed to create, deploy, and execute various application programs on various platforms by packaging the necessary dependencies of the application [18].

Figure 1 shows the overview of the UPC system. The UPC system adopts the *masterworker model*. Users submit computing jobs to the *UPC master* through the *UPC web server*. After synchronizations of the jobs, the *UPC master* assigns the submitted jobs to the appropriate *UPC workers*. Each worker computes its assigned jobs and returns the results to the master upon completion. Users can access the results at the web browser.

**Figure 1.** Overview of UPC system.

For further details, the usage flow of the UPC system will be described:


#### *3.2. OpenPose*

Next, we review *OpenPose*. It has been developed by researchers at Carnegie Mellon University and is an popular open-source software for real-time human pose estimation [3]. It extracts the feature points, called *keypoints*, of the human body in the given image using *Convolutional Neural Network (CNN)*. The *keypoints* represent the important joints in a human body, the contours of eyes, lips in the face, fingertips, and joints in the hands and feet. Using the *keypoints*, the shapes of a body, face, hands, and feet can be described. Since it has been developed based on CNN, the CPU time is very long when computed on a conventional PC.

*OpenPose* is used in our group for developing the *exercise and performance learning assistant system (EPLAS)* to assist practicing exercises or learning performances by themselves at home [19]. *EPLAS* offers video content of *Yoga* poses by instructors whose performances should be followed by users. During the practice, it automatically takes photos of important scenes of the user. Then, it extracts the keypoints of the human body using *OpenPose* to rate the poses in the photos by comparing the coordinates of them between the user and the instructor.

#### *3.3. OpenFOAM*

Then, we review *OpenFOAM*. It is an open-source software for the *computational fluid dynamics (CFD)* simulations and has been developed primarily by *OpenCFD Ltd.* (Bracknell, UK) It has an extensive range of features to solve anything from complex fluid flows involving chemical reactions, turbulence, and heat transfer, to acoustics, solid mechanics, and electromagnetics [4]. Furthermore, the optimal parameter selection is critical for the high accuracy of the results, and it needs a lot of iterations of selecting parameters in *OpenFOAM* and running it with the parameter values. We applied the *parameter optimization method* for *OpenFOAM* [20]; it needs to run *OpenFOAM* with a lot of different parameters.

Meanwhile, it is also applied for developing the *air conditioning guidance system* [21] in our research. The estimation or prediction of the distributions of the temperature or humidity inside a room using this simulation model is necessary to properly control the air conditioner. By estimating the room environment changes under various actions, it will be possible to decide when the air conditioner is turned on or off. Even the timing to open or close windows in the room can be selected. To estimate or predict the distributions in a room together with sensors, the CFD simulation using *OpenFOAM* has been investigated. Then, the optimization of the parameters in *OpenFOAM* is critical in order to fit the simulation results well with the corresponding measured ones.

#### *3.4. Code Testing*

Finally, we review the *code testing* in the *Android programming learning assistance system (APLAS)*. *APLAS* has been developed in our group as the automatic and self-learning system for Android programming using *Java* and *XML* [5,6]. The *code testing* is the process to validate a source code by running the corresponding *test code* on a testing framework. To confirm the validity of the answer source code from a student in satisfying the required specifications in the assignment, *APLAS* implements the *code testing* function using *JUnit* for unit testing of *Java* codes [22] and *Robolectric* for integration testing with *XML* codes [23,24]. *APLAS* needs to run the *code testing* function with a lot of different source codes from many students, which usually takes a long time.

In *ALPAS*, *Java* codes can be directly tested on *JUnit*. However, the Android-specific components, such as the *Layout*, the *Activity*, the *Event Listener*, and the *Project Resources* that will be described in *XML*, cannot be directly tested on *JUnit*. The building tool *Gradle* is used to build and integrate them as *Java* classes. Then, *Robolectric* is used to generate *Java* objects—called *shadow objects*—for them, so that they can be tested on *JUnit*.

#### **4. Proposal of Static Uniform Job Assignment Algorithm**

In this section, we present the static *uniform job* assignment algorithm to workers in the UPC system.

#### *4.1. Objective*

To design the algorithm, it is observed that when the *makespan* of every worker becomes equal, the objective of the problem on the *makespan* minimization can be achieved. Otherwise, the maximum *makespan* can be reduced by moving some jobs at the bottleneck worker which determines this maximum *makespan* to other workers, if the number of assigned jobs to any worker can take a real number. Only when every worker has the same *makespan*, the maximum *makespan* cannot be reduced.

$$\min \text{imizar}\{\max(m\_w^t)\} \text{ for } t \in T, w \in \mathcal{W} \tag{1}$$

The minimization of the maximum *makespan* among all the workers is given as the objective of the problem, where *makespan m<sup>t</sup> <sup>w</sup>* at worker *w* for type *t* is given by the summation of the CPU time for preparation and execution.

#### *4.2. Simultaneous Linear Equations*

In this paper, the following *simultaneous linear equations* have been derived to find the optimal job-worker assignment, such that the estimated CPU time required to complete the assigned jobs becomes equal among all the workers. The solutions of the *simultaneous linear equations* will be the *lower bound* on *makespan*. Since the solutions become real numbers in general, the integer number of assigned jobs to each worker should be introduced to them.

$$\begin{aligned} \mathbf{C}\_i^t + \frac{\mathbf{R}\_{i,D\_i}^t}{D\_i} \times \mathbf{x}\_i^t &= \mathbf{C}\_j^t + \frac{R\_{j,D\_j}^t}{D\_j} \times \mathbf{x}\_j^t \\ \text{for } i \neq j, i \in \mathcal{W}, j \in \mathcal{W}, t \in T. \end{aligned} \tag{2}$$

To satisfy the objective of the equal CPU time among the workers, *R<sup>t</sup> <sup>w</sup>*,*Dw* /*Dw* gives the best CPU time to solve one job at worker *w* by running *Dw* jobs.

#### *4.3. Problem Formulation*

To present the static *uniform jobs* assignment algorithm to workers in the UPC system, the problem to be solved is formulated here.

#### 4.3.1. Variables

The following variables are defined for the problem to be solved:


#### 4.3.2. Constants

The following constants are given as the inputs to this problem:


Here, *Dw* represents the number of simultaneously running jobs for job type *t* at worker *w*, which maximizes the number of completed jobs per unit time. This is constant for any job type in each application, because it depends on the common program in the application for every job type.

*Ct <sup>w</sup>* represents the CPU time required to initiate the execution of the program at worker *w*. For example, in the *code testing* application, it represents the CPU time to initiate the *Gradle Wrapper* daemon and generate *shadow objects* that are necessary to run the *code testing function*.

*Rt <sup>w</sup>*,*<sup>d</sup>* can be measured using any worker by running jobs for job type *t* while increasing the number of running jobs in parallel from 1 until *Dw*.

#### 4.3.3. Constraints

The following two constraints must be satisfied in the problem:

• The total number of the assigned jobs to workers must be equal to *N<sup>t</sup>* for any type *t*.

$$\sum\_{w \in \mathcal{W}} \mathbf{x}\_w^t = \mathbf{N}^t \ (t \in T) \tag{3}$$

• Any worker cannot run *d* jobs in parallel when *d* is larger than the *Dw* (let *dw* for worker *w*) due to the PC specifications.

$$d\_{w} \le D\_{w} \tag{4}$$

#### *4.4. Conditions for Uniform Job Assignment*

For the *uniform job* assignment to workers in the UPC system, the following conditions are assumed:


#### *4.5. Static Uniform Job Assignment Algorithm*

Here, we note that the CPU time may be different depending on the number of running jobs in parallel in each worker that has multiple cores. To reduce the CPU time by increasing the job completion throughput, *Dw* jobs of type *t* should run at worker *w* as much as possible, since it will give the best throughput. Based on this observation, we present the three-step static *uniform job* assignment algorithm. Figure 2 shows the flowchart of the proposal.

#### 4.5.1. First Step

By solving the *simultaneous linear equations* composed of (2) and (3), the optimal number of assigned jobs of type *t* to worker *w*, ˆ*x<sup>t</sup> <sup>w</sup>*, is obtained, assuming that any real value is acceptable for it.

**Figure 2.** Flowchart of the proposal.

#### 4.5.2. Second Step

The solution in the first step becomes feasible only when ˆ*x<sup>t</sup> <sup>w</sup>* is a multiple of *Dw* for type *t*. Unfortunately, ˆ*x<sup>t</sup> <sup>w</sup>* does not satisfy the condition, in general. Therefore, in the second step, as the closest integer number to satisfy the condition, the following ˜*x<sup>t</sup> <sup>w</sup>* jobs will be assigned to the worker (worker *w*), where *y* gives the largest integer equal to or smaller than *y*:

$$\hat{\mathfrak{x}}\_w^\dagger = \lfloor \frac{\mathfrak{x}\_w^\dagger}{D\_w} \rfloor \times D\_w \tag{5}$$

Then, the number of the remaining jobs (let *r<sup>t</sup>* for type *t*) is calculated by:

$$r^t = N^t - \sum\_{w \in \mathcal{W}} x\_w^{\overline{t}} \tag{6}$$

Besides, the estimated *makespan* for each worker (let *em<sup>t</sup> <sup>w</sup>* for worker *w* and job type *t*) after the job assignment is calculated by:

$$em\_w^t = \mathbb{C}\_w^t + R\_{w, D\_w}^t \times \frac{\bar{\mathfrak{X}}\_w^t}{D\_w} \tag{7}$$

Therefore, after completing the procedures for all the job types, the estimated *makespan* for each worker is calculated by:

$$EM\_w = \sum\_{t \in T} em\_w^t \tag{8}$$

As the objective of the algorithm, the maximum estimated *makespan* among the workers is calculated by:

$$EM = \{ \max(EM\_w) \} \text{ for } w \in \mathcal{W} \tag{9}$$

#### 4.5.3. Third Step

In the third step, the remaining jobs (*rt*) in the second step will be assigned to workers in a greedy way, such that the increase in the maximum estimated *makespan EM* is minimized. It is noted that the remaining jobs may exist for any job type. Here, to utilize the parallel job computation using multiple threads on multiple cores for each worker as much as possible, the simultaneous assignment of multiple jobs to one worker is always considered.

1. Find the worker whose *EM*ˆ *<sup>w</sup>* is smallest among the workers (let worker *w*).

$$E\hat{\mathcal{M}}\_{w} = E\mathcal{M}\_{w} + \mathcal{R}\_{w,D\_{w}}^{t} \tag{10}$$

2. Assign Δ*x<sup>t</sup> <sup>w</sup>* jobs to worker *w*.

$$
\Delta \mathbf{x}\_w^t = \begin{cases} D\_{w\prime} & r\_t > D\_w \\ r\_{t\prime} & r\_t \le D\_w \end{cases} \tag{11}
$$

3. Update the number of the remaining jobs (*rt*), and the number of assigned jobs and *makespan* of the worker *w* by:

$$\begin{aligned} \mathbf{x}\_{\textit{w}}^{t} &= \mathbf{x}\_{\textit{w}}^{t} + \Delta \mathbf{x}\_{\textit{w}\prime}^{t} \\ \mathbf{E}M\_{\textit{w}} &= EM\_{\textit{w}} + R\_{\textit{w},\Delta \mathbf{x}\_{\textit{w}}^{t}\prime}^{t} \\ r\_{t} &= r\_{t} - \Delta \mathbf{x}\_{\textit{w}}^{t} \end{aligned} \tag{12}$$

4. If the number of the remaining jobs becomes zero (*rt* = 0), terminate the procedure.

5. Go to 1.

#### **5. Evaluation**

In this section, we evaluate the proposal through extensive experiments which are running jobs in three applications on the testbed UPC system.

#### *5.1. Testbed UPC System*

Table 1 shows the PC specifications in the testbed UPC system. One master and six workers are used here.

**Table 1.** PC specifications.


#### *5.2. Jobs*

Table 2 shows the specifications of the jobs for the eight job types in our experiments. For the *code testing* application in *APLAS*, six job types are prepared, where each job type represents one assignment to students in *APLAS*. These job types run the same programs of *JUnit* and *Robolectric*, but accept many different data of answer source codes and test codes. For the other applications, only one job type is considered.



#### *5.3. CPU Time*

Table 3 shows the constant CPU time required to start running the jobs on each worker for each of the six job types. Tables 4–6 show the increasing CPU time when the number of jobs is increased by one until the number for the best throughput for each type.

Through preliminary experiments, we found the number of simultaneously running jobs for the highest throughput for each worker. For *code testing* in *APLAS*, *PC1*, *PC2*, and *PC3* can run only one job in parallel due to the low specifications. This number is two for *PC4*, five for *PC5*, and six for *PC6*. For *OpenPose*, any worker can only execute one job because it uses a lot of threads to compute CNN. For *OpenFOAM*, for each worker, the CPU time is constant at any number of simultaneously running jobs until it reaches the number of cores in the worker.

**Table 3.** Constant CPU time to start jobs (s).


#### **Table 4.** Increasing CPU time at *PC1*∼*PC4* (s).


**Table 5.** Increasing CPU time at *PC5* (s).


**Table 6.** Increasing CPU time at *PC6* (s).


#### *5.4. Comparative Algorithms*

For performance comparisons, we implemented two simple algorithms to assign *non-uniform jobs* to workers.

The first one is the *First-Come-First-Serve (FCFS)* algorithm. It assigns each job to the first available worker, starting from the worker with the highest specification until the one with the lowest. It limits the worker to executing only one job at a time.

The second is the *best throughput-based FCFS (T-FCFS)* algorithm. The difference between T-FCFS and *FCFS* is that each worker may execute multiple jobs simultaneously until the best throughput.

#### *5.5. Total Makespan Results*

Table 7 compares the maximum *makespan* results for each job type when the testbed UPC system runs the jobs by following the assignments found by the algorithms. Furthermore, it shows the *lower bound (LB)* on the maximum *makespan* found at *First Step* of the proposed algorithm for the reference of them.


**Table 7.** Maximum *makespan* results (s).

The results indicate that for any job type, the maximum *makespan* result by the proposal is better than the results by the two compared algorithms and is close to the lower bound. Thus, the effectiveness of the proposal is confirmed. It is noted that the results by *FCFS* are far larger than the ones by the others because FCFS does not consider simultaneous multiple job executions for a worker.

#### *5.6. Individual Makespan Results*

For reference, Tables 8–10 show *makespan* or the total CPU time of each worker and the largest CPU time difference between the workers and the three algorithms. For *OpenFOAM*, no job was assigned to *PC1*–*PC4*, because all of the 32 jobs can be executed simultaneously at *PC5* and *PC6*. The largest CPU time difference by the proposal is smaller than the ones by the others, except for *ColorGame*, *SoccerMatch*, *AnimalTour*, and *MyLibrary*, where in Table 4, the increasing CPU time of *PC1* is much larger than other workers, and the far smaller number of jobs was assigned. Therefore, the proposal can balance well the job assignments among the workers.




**Table 9.** T-FCFS *makespan* detail (s).

**Table 10.** Proposal *makespan* detail (s).


#### *5.7. Discussions*

The results in Table 7 show improvements of maximum *makespan* results by the proposed algorithm if compared with *T-FCFS*. However, some differences can be observed against the *lower bound*.

The current algorithm can find the assignment of some remaining jobs to workers, and assign an integer number of jobs to any worker in a greedy way, after the real number solutions are obtained by solving the *simultaneous linear equations*. A greedy method is usually difficult to give a near-optimum solution, since it only considers the local optimality under the current assignment.

To improve the solution quality, a local search method using iterations has often been adopted for solving combinatorial optimization problems, including this study. Therefore, we will study the use of a local search method for the remaining job assignment in the proposed algorithm.

#### **6. Extension to Multiple Job Types Assignment**

In this section, we extend the proposed algorithm to the case when jobs for multiple job types are assigned together.

#### *6.1. Algorithm Extension*

In *First Step* of the proposed algorithm, the linear equations are modified in this extension to consider the CPU time to complete all the jobs for the plural job types assigned to each worker:

$$\begin{aligned} \sum\_{t \in T} \left( \mathbf{C}\_i^t + \frac{R\_{i, D\_i}^t}{D\_i} \times \mathbf{x}\_i^t \right) &= \sum\_{t \in T} \left( \mathbf{C}\_j^t + \frac{R\_{j, D\_j}^t}{D\_j} \times \mathbf{x}\_j^t \right) \\ \text{for } i \neq j, i \in \mathcal{W}, j \in \mathcal{W}. \end{aligned} \tag{13}$$

The number of variables to be solved is |*W*||*T*|, where |*W*| represents the number of workers and |*T*| represents the number of job types, respectively. Thus, |*W*||*T*| linear equations are necessary to solve them. In the original algorithm, for each job type, (|*W*| − 1) linear equations are derived for the CPU time equality and one equation is for the job number. Thus, |*W*||*T*| equations can be introduced.

However, in this extension, the total number of linear equations for the CPU time equality is reduced to (|*W*| −1) because all the job types need to be considered together here. Therefore, to solve the linear equations uniquely, the following (|*W*| − 1)(|*T*| − 1) linear equations will be introduced by considering the total CPU time for (|*T*| − 1) job types together for (|*T*| − 1) combinations of (|*T*| − 1) job types, in addition to the total CPU time for |*T*| job types together in (13):

$$\sum\_{t \in T - \{u\}} (\mathbf{C}\_i^t + \frac{R\_{i, D\_i}^t}{D\_i} \times \mathbf{x}\_i^t) = \sum\_{t \in T - \{u\}} (\mathbf{C}\_j^t + \frac{R\_{j, D\_j}^t}{D\_j} \times \mathbf{x}\_j^t) \tag{14}$$
  $for \ i \neq j, i \in \mathcal{W}, j \in \mathcal{W}, u \in T.$ 

where *T* − {*u*} represents the set of the job types in *T* except for job type *u*.

The (|*T*| − 1) combinations of (|*T*| − 1) job types are selected by excluding the combination where the following estimated total CPU time to execute all the jobs in the remaining job types on *PC6* is smallest:

$$\sum\_{t \in T - \{u\}} \left( \mathbf{C}\_6^t + \frac{R\_{6, D\_6}^t}{D\_6} \times N^t \right) \tag{15}$$

Then, in *Second Step* and *Third Step*, the estimated *makespan* for each worker and the maximum estimated *makespan* among the workers are modified to consider all the given job types together.

#### *6.2. Total Makespan Results*

Table 11 shows the maximum *makespan* results when the testbed UPC system runs the jobs by following the assignments by the extended algorithm. When compared with the result by the original algorithm, it is reduced by 5%, and becomes closer to the *lower bound*. The difference between our result and the lower bound is very small. Thus, this extension is effective when plural job types are requested at the UPC system together.

**Table 11.** Maximum *makespan* results (s) by proposal.


#### *6.3. Discussions*

The result in Table 11 confirms some reduction in the total *makespan* result by the extended algorithm. However, there is still a difference when compared to the *lower bound*. Thus, it is necessary to further improve the algorithm.

One idea for this improvement in the extended algorithm will not be to limit the exclusion of one job type combination—where the estimated total CPU time to execute all jobs in the remaining job types on *PC6* is the smallest—and to generate the linear equations for the CPU time equality. Instead, every combination will be excluded one by one to obtain the result for each combination exclusion. Then, the best one will be selected among them.

#### **7. Conclusions**

This paper proposed the static *uniform job* assignment algorithm to workers in the UPC system. The *simultaneous linear equations* have been derived to find the optimal assignment of minimizing the maximum *makespan* among the workers, where the CPU time to complete the assigned jobs becomes equal among all the workers.

For an evaluation, the 651 *uniform jobs* in three applications, *OpenPose*, *OpenFOAM*, and *code testing* in *APLAS*, were considered to run on six workers in the testbed UPC system, and the *makespan* was compared with the results by two simple algorithms and the lower bounds. The comparisons confirmed the effectiveness of the proposal.

The novelty of the proposal is that with a very simple formula, it is able to provide the near-optimal solutions to *NP-complete* problems in the *User-PC computing (UPC) system*, a typical distributed system. The current algorithm limits the jobs whereby the computing time for a worker is nearly equal. This limitation can simplify our approach of considering the simple assignment of the number of jobs for each worker without considering the differences among individual jobs.

Fortunately, it is possible to alleviate this limitation by considering the granularity of the CPU time for a worker. The CPU time of a job in suitable applications to the proposal is often proportional to the number of iteration steps before the termination or the number of elements in the model. By considering a multiple of a constant number of iteration steps, the CPU time can be estimated even if the number of iteration steps is widely changed with this granularity; this finding will be in future studies.

In future studies, we will also improve the algorithm for remaining job assignments and simultaneous job assignments of multiple job types, and we will study the combination of *uniform jobs* and *non-uniform jobs* in the job-worker assignment algorithm for the UPC system.

**Author Contributions:** Conceptualization, N.F.; Data curation, X.Z. and H.H.; Resources, H.H., A.K., I.T.A., Y.H. and Y.W.S.; Software, X.Z.; Supervision, N.F.; Writing–original draft, X.Z.; Writing–review & editing, N.F. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Data Availability Statement:** Not applicable.

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

#### **Abbreviations**

The following abbreviations are used in this paper:


#### **References**


### *Article* **The Assignment Problem and Its Relation to Logistics Problems**

**Milos Seda**

Institute of Automation and Computer Science, Faculty of Mechanical Engineering, Brno University of Technology, Technicka 2896/2, 623 00 Brno, Czech Republic; seda@fme.vutbr.cz

**Abstract:** The assignment problem is a problem that takes many forms in optimization and graph theory, and by changing some of the constraints or interpreting them differently and adding other constraints, it can be converted to routing, distribution, and scheduling problems. Showing such correlations is one of the aims of this paper. For some of the derived problems having exponential time complexity, the question arises of their solvability for larger instances. Instead of the traditional approach based on the use of approximate or stochastic heuristic methods, we focus here on the direct use of mixed integer programming models in the GAMS environment, which is now capable of solving instances much larger than in the past and does not require complex parameter settings or statistical evaluation of the results as in the case of stochastic heuristics because the computational core of software tools, nested in GAMS, is deterministic in nature. The source codes presented may be an aid because this tool is not yet as well known as the MATLAB Optimisation Toolbox. Benchmarks of the permutation flow shop scheduling problem with the informally derived MIP model and the traveling salesman problem are used to present the limits of the software's applicability.

**Keywords:** assignment problem; traveling salesman problem; vehicle routing problem; flow shop scheduling problem; GAMS, genetic algorithm

#### **1. Introduction**

The *Assignment Problem* (abbreviated to AP) [1] and its mathematical model is a problem that is the basis of the field of combinatorial optimization [2,3]. The problem in its basic form has been successfully handled by the discovery of Harold Kuhn, who proved that his method [4], derived from the results of the theoretical work of the Hungarian mathematicians Dénes König and Jenö Egerváry and dubbed the *Hungarian method* in their honor, finds a solution in polynomial time O(*n*3) [5] in an efficient implementation.

However, this does not make the assignment problem less interesting because it has many analogs in bipartite graph matching problems of the graph theory [5,6].

The assignment problem is most extensively addressed in [5,7], where we find theoretical foundations for the existence of perfect matching, implementation details for the Hungarian method, and a number of other related problems such as the *k*-cardinality assignment problem, the semi-assignment problem, the bottleneck assignment problem, the algebraic assignment problem, quadratic assignment problems, and multi-index assignment problems.

These are still variants closely related in meaning to the basic version of the matching problem, although the time complexity may no longer be polynomial, and, in the case of the quadratic matching problem, it is a non-linear problem.

In this paper, however, we focus on similarities of a different kind, namely in the mathematical model, where a completely different relationship to the underlying problem may be found because it allows, often with small modifications, for expressing problems from a different application domain, to switch between linear, mixed integer, and even non-linear problem classes, thus changing its computational complexity and solvability.

Perhaps the most well-known combinatorial optimization problem, the Travelling Salesman Problem [8], is a slight modification of this, with the Vehicle Routing Problem [9]

**Citation:** Seda, M. The Assignment Problem and Its Relation to Logistics Problems. *Algorithms* **2022**, *15*, 377. https://doi.org/10.3390/a15100377

Academic Editor: Frank Werner

Received: 7 September 2022 Accepted: 13 October 2022 Published: 16 October 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

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

also stemming from this. On the other hand, from the assignment problem, with other modifications, we can easily move to logistic distribution operations [10], agricultural applications [11], set covering problems [12] with many interesting applications in telecommunications [13], and scheduling problems [14].

Concerning two selected problems of exponential complexity, the Travelling Salesman Problem (TSP) and the Permutation Flow Shop Scheduling Problem (PFSSP) [15], we will also deal with their solvability. We have very good experience in solving sets covering problems of O(2*n*) complexity [13] using GAMS (General Algebraic Modelling System), where it is possible to find the optimal solution in the available time even for instances with matrices of hundreds of rows and thousands of columns, and it has also proven itself in solving the problem of finding the Steiner minimum tree in networks, which also has exponential time complexity.

Since the two problems mentioned above are permutational in nature with the factorial time complexity, they are more challenging than the set covering problems.

For extremely large TSP instances (many hundreds of cities), heuristics must be used, e.g., differential evolution [16], genetic algorithm [17–20], memetic search [21], simulated annealing [22], neural network [23], and improved neighbourhood search algorithms [24,25].

Many stochastic heuristics are inspired by the behaviour of animals in nature, e.g., deer [26], spider monkey [27], hyena [28], wolf [29], cuckoo [30,31], sparrow [32], frog [33], and ant colony [34].

On the other hand, stochastic heuristic methods are not suitable for TSP instances up to 100 cities because they may not find the optimal solution and the convergence time is often unsure, as, e.g., shown in the comparison of different methods in [28].

However, there are also approaches based on deterministic methods such as cutting plane [35], branch and bound [36] and branch and cut [37,38].

According to listings in the GAMS environment, the latter method is in some way incorporated into GAMS and, therefore, it makes sense to explore its limits of applicability. These, together with the GAMS source code, are discussed in detail in Section 6.

Scheduling problems seem far from the assignment problem. But one of them, the PF-SSP, shares with the assignment problem a permutational nature in the ordering of jobs, where each job (with its operations) is assigned to exactly one position and each position can contain only one job (with its operations), which corresponds in the assignment problem to the fact that each task is assigned to a single worker and each worker solves only one task. There are additional constraints, and the aim is to minimize the total scheduling time (makespan), but we can still say that the derived PFSSP model is related to the assignment problem.

As in the case of the TSP, heuristic methods are used for large instances of different variants of flow shop scheduling problems, e.g., differential evolution [39], genetic algorithm [40,41], genetic programming [42], memetic algorithm [43], tabu-search [44], harmony search [45], iterated greedy algorithms [46–48], multi-local search [49], hybrid metaheuristics [50,51], reinforcement learning [52], fireworks algorithm [53], and also nature-inspired algorithms, e.g., ant colony optimization [54], firefly particle swarm optimization [55], migrating birds optimization [56] and whale swarm algorithm [57].

However, the exact methods [58,59], linear programming approach [60] and branch and bound [61], are also applicable so we will again focus on the usability of the GAMS tool.

#### **2. The Assignment Problem Model**

In most common problem formulation, we have *n* workers who need to be assigned *n* tasks in such a way that each worker is assigned a single task and each task is solved by a single worker.

For each worker-task pair, we know the time it takes the worker to complete the task. The task is to find an assignment that minimizes the total time to complete all tasks.

Let *cij* denote the time taken by the *i*th worker for the *j*th task. The decision variables are binary, *xij* = 1 if the *i*th worker is assigned the *j*th task, *xij* = 0 in the opposite case. Then, the problem can be formulated as follows:

$$z = \sum\_{i=1}^{n} \sum\_{j=1}^{n} c\_{ij} x\_{ij} \to \min \tag{1}$$

subject to

$$\sum\_{i=1}^{n} \mathbf{x}\_{ij} = 1, \ j = 1, \dots, n$$

$$\sum\_{j=1}^{n} x\_{ij} = 1, \ i = 1, \ \dots, n$$

$$\mathbf{x}\_{i\bar{j}} \in \{0, 1\}, \ i = 1, \ldots, n, \ j = 1, \ldots, n \tag{4}$$

Equation (2) ensures that each task is assigned to a single worker, and Equation (3) ensures that each worker is assigned a single task.

The assignment problem can also be viewed as a problem of finding a permutation

$$
\begin{pmatrix} 1 & 2 & \dots & n \\ \pi\_1 & \pi\_2 & \dots & \pi\_n \end{pmatrix},
$$

where *i*th worker is assigned to task *π<sup>i</sup>* and

$$z = \sum\_{i=1}^{n} c\_{i\pi\_i} \to \min$$

Since the number of different permutations of *n* elements is *n*!, it is not possible to find the optimal solution for large instances in the available time by enumerating all possibilities. However, due to the Hungarian method mentioned above, we no longer use this approach.

#### **3. Routing Problems**

With a different interpretation of the variables and a possible extension of the constraints, the assignment problem changes into a series of other problems. In this section, we consider two routing problems.

#### *3.1. Travelling Salesman Problem*

The *Travelling Salesman Problem* (TSP) [8,62] is mathematically similar to the assignment problem model, differing only in one additional constraint, but the meaning of the decision variables *xij* = 0 is different. It is formulated as follows: Given *n* cities and distances among them, the objective is to find a round trip through all cities with a minimum length (alternatively, with a minimum total transportation cost).

Since the starting city 1 is fixed, the number of routes is given by the permutations of cities 2, 3, ..., *n*, and is therefore equal to (*n* − 1)!. If there are no one-way segments anywhere in the transportation between cities, routes in reverse order of cities do not affect the length, and then we can reduce the number of routes to (*n* − 1)!/2, but still the time complexity of exploring all routes is O(*n*!).

If we denote by *cij* the distance between cities *i* and *j* (alternatively, the price of transportation between cities *i* and *j*), *xij* a binary decision variable that takes the value 1 when city *j* on the route immediately follows city *i*, otherwise it takes the value 0, *δ<sup>i</sup>* is the order of city *i* on the route, then the Travelling Salesman Problem can be formulated as follows:

$$z = \sum\_{i=1}^{n} \sum\_{j=1}^{n} c\_{ij} \mathbf{x}\_{ij} \to \min \tag{5}$$

subject to

$$\sum\_{i=1}^{n} \mathbf{x}\_{ij} = 1, \ j = 1, \dots, n$$

$$\sum\_{j=1}^{n} x\_{ij} = 1, \ i = 1, \ldots, n \tag{7}$$

$$
\delta\_i - \delta\_j + n\mathbf{x}\_{ij} \le n - 1, \text{ i } \ne j, \text{ i } = 2, \dots, n, \text{ j } = 2, \dots, n \tag{8}
$$

$$\mathbf{x}\_{i\mathbf{j}} \in \{0, 1\}, \ i = 1, \ldots, n, \ j = 1, \ldots, n \tag{9}$$

Constraints (6) and (7) ensure that each city (vertex of the graph) is traversed exactly once (entered and left exactly once); the system of subtour elimination constraints (8), referred to in the English literature by Miller, Tucker, and Zemlin as *MTZ constraints*, prevents the formation of subtours, as we will show below in Theorem 1.

Equation (8) follows from the following reasoning:

(i) If *xij* = 1, then *j* is the immediate successor of *i*, and if *δ<sup>i</sup>* = *t*, then *δ<sup>j</sup>* = *t* + 1. Hence, *δ<sup>i</sup>* − *δ<sup>j</sup>* + *nxij* = *t* − (*t* + 1) + *n* = *n* − 1.

(ii) If *xij* = 0, then *δ<sup>i</sup>* − *δ<sup>j</sup>* + *nxij* = *δ<sup>i</sup>* − *δ<sup>j</sup>* and this difference in the order of the cities in the route for *i* = *j* can be at most equal to *n* − 1.

From (i) and (ii), a common conclusion *δ<sup>i</sup>* − *δ<sup>j</sup>* + *nxij* ≤ *n* − 1 already follows, which, for all combinations of feasible values of *i* and *j*, is expressed by inequality (8).

Without the constraint (8), constraints (6) and (7) are satisfied by splitting the route into several subtours, e.g., for 15 vertices, the two conditions mentioned above are satisfied by the subtours 1 − 3 − 7 − 9 − 12 − 1, 2 − 4 − 10 − 11 − 13 − 15 − 2 and 5 − 6 − 8 − 14 − 5.

**Theorem 1.** (Miller, Tucker, Zemlin) *The variables xij* ∈ {0, 1}, *i* = 1, ... , *n*, *j* = 1, ... , *n satisfying constraints* (6) *and* (7) *form a Hamiltonian circle if and only if the subtour elimination constraints* (8) *are satisfied.*

**Proof.** Suppose that *xij* satisfies the subtour elimination constraints but does not form a Hamiltonian circle. Then *xij* due to (6) and (7) form at least two subtours, one containing the initial vertex 1 and another without it. Let *S* be a subtour that does not contain vertex 1 and let *E*(*S*) be the set of edges in *S*. Summing the conditions over the edges of *E*(*S*) we get:

$$\sum\_{(i,j)\in E(S)} (\delta\_i - \delta\_j + n\mathfrak{x}\_{ij}) \le (n-1)|E(S)|\_{\prime}$$

since the values of *δ<sup>i</sup>* and *δ<sup>j</sup>* eliminate each other in this subtour, we get

$$n|E(S)| \le (n-1)|E(S)|,$$

which is a contradiction.

Assume now that *xij* forms a Hamiltonian circle. If 1 is the initial vertex of this circle, and for each vertex *i* = 1, *δ<sup>i</sup>* = *k*, if *i* is the *k*th vertex of the Hamiltonian circle, then it is clear that the conditions (8) are satisfied.

#### *3.2. Vehicle Routing Problem*

A generalization of the Travelling Salesman Problem is the *Vehicle Routing Problem* (VRP) [9,63–66] where a desired quantity of goods needs to be delivered from a central depot to customers by vehicles of a certain capacity.

We are looking for closed routes of individual vehicles that start and end at the depot, each customer is served exactly once by exactly one vehicle, the requirements of all customers are met and the total transport costs are minimal.

Consider the following notation:

*n* . . . number of customers

0 . . . depot (start and end of each vehicle's route)

*K* . . . number of (identical) vehicles


*xij* ... binary decision variable equal to 1 if *j* is immediately followed by *i* on the route, *xij* = 0 otherwise

*δ<sup>i</sup>* . . . the load left in the vehicle after visiting customer *i*.

$$z = \sum\_{i=0}^{n} \sum\_{j=0}^{n} c\_{ij} \mathbf{x}\_{ij} \to \min \tag{10}$$

subject to

$$\sum\_{i=0}^{n} \mathbf{x}\_{i\bar{j}} = 1, \ j = 1, \ldots, n \tag{11}$$

$$\sum\_{j=0}^{n} x\_{ij} = 1, \ i = 1, \ \dots, n \tag{12}$$

$$\sum\_{i=1}^{n} \mathbf{x}\_{i0} = K \tag{13}$$

$$\sum\_{j=1}^{n} \varkappa\_{0j} = K \tag{14}$$

$$0 \le \delta\_i \le Q - d\_i, \ i = 1, \dots, n \tag{15}$$

$$
\delta\_i - \delta\_j + Q \mathbf{x}\_{ij} \le Q - d\_j, \text{ i } \ne j, \text{ i } = 1, \dots, n, \ j = 1, \dots, n, \text{ such that } d\_i + d\_j \le Q \tag{16}
$$

$$\text{tr}\_{ij} \in \{0, 1\}, \ i = 0, \ldots, n, \ j = 0, \ldots, n \tag{17}$$

In the model, (11) and (12) ensure that exactly one vehicle arrives at each customer (11) and exactly one vehicle leaves it (12). Equations (13) and (14) ensure that all *K* vehicles return to the depot (13) and all *K* vehicles leave the depot (14).

Equation (16) is analogous to the MTZ constraints in the Travelling Salesman Problem preventing the formation of partial circuits and at the same time ensuring that the requirements of customers *i* and *j* can be met when traveling from *i* to *j* [67].

The Vehicle Routing Problem has many other specific formulations, e.g., there may be a larger number of depots available, and customers are only ready to receive delivery of goods at certain time intervals. For more details, see the sources listed at the beginning of this section.

#### **4. Distribution Problems**

Distribution problems have many different formulations, first, we consider the classical Hitchcock's *Transportation/Transshipment Problem* with *m* suppliers (sources, warehouses) and *n* customers (consumers), where we assume the transportation of a single type of material (goods) with an objective to minimize the total cost of transporting the material [68].

Assume the following notation:

*ai*, *i* = 1, . . . , *m* ... capacity (stocks) of suppliers,

*bj*, *j* = 1, . . . , *n* ... customer requirements,

*cij*, *i* = 1, ... , *m*, *j* = 1, ... , *n* ... the matrix of rates for the transport of a unit quantity between the *i*th supplier and the *j*th customer,

*xij*, *i* = 1, ... , *m*, *j* = 1, ... , *n* ...the sought quantity transported between the *i*th supplier and the *j*th customer.

If total stocks are equal to total requirements, this means:

$$\sum\_{i=1}^{m} a\_i = \sum\_{j=1}^{n} b\_{j\prime} \tag{18}$$

we are talking about a *balanced distribution problem*, where all stocks are exhausted and all demands are met, and the following mathematical model corresponds to this:

$$z = \sum\_{i=1}^{m} \sum\_{j=1}^{n} c\_{ij} \mathbf{x}\_{ij} \to \min \tag{19}$$

subject to

$$\sum\_{j=1}^{n} x\_{ij} = a\_i, \ i = 1, \ \dots, \ m \tag{20}$$

$$\sum\_{i=1}^{m} x\_{ij} = b\_{j\prime} \; j = 1, \; \dots \; \; n \tag{21}$$

$$\mathbf{x}\_{i\bar{j}} \ge \mathbf{0}, \; i = 1, \; \dots, m, \; j = 1, \; \dots, n \tag{22}$$

Equation (20) corresponds to the stock drawdown, and Equation (21) expresses the fulfillment of requirements.

Obviously, the assignment problem is a special case of the balanced transportation problem, where:

$$m = n$$

$$a\_i = 1, \ i = 1, \ \dots, m$$

$$b\_j = 1, \ j = 1, \ \dots, n$$

$$x\_{ij} \in \{0, 1\}, \ i = 1, \ \dots, m, \ j = 1, \ \dots, n$$

However, the balanced case is rare in practice, usually, total stocks exceed total requirements, i.e., *<sup>m</sup>*

$$\sum\_{i=1}^{m} a\_i > \sum\_{j=1}^{n} b\_j \tag{23}$$

In this case, all requirements can be met, but not every stock will be used up. The model then changes as follows:

$$z = \sum\_{i=1}^{m} \sum\_{j=1}^{n} c\_{ij} x\_{ij} \to \min \tag{24}$$

subject to

$$\sum\_{j=1}^{n} x\_{ij} \le a\_{i\prime} \text{ i } i = 1 \text{ } \dots \text{ } m \tag{25}$$

$$\sum\_{i=1}^{m} x\_{ij} = b\_{j\prime} \; j = 1, \; \dots, n \tag{26}$$

$$\mathbf{x}\_{ij} \ge \mathbf{0}, \ i = 1, \ldots, m, \ j = 1, \ldots, n \tag{27}$$

In the case of material shortages, the opposite situation may occur, where the total stock is insufficient for the total requirements, i.e.,

$$\sum\_{i=1}^{m} a\_i < \sum\_{j=1}^{n} b\_j \tag{28}$$

This means that stocks are used up but not all requirements can be met. The model must then be modified as follows:

$$z = \sum\_{i=1}^{m} \sum\_{j=1}^{n} c\_{ij} \mathbf{x}\_{ij} \to \min \tag{29}$$

subject to

$$\sum\_{j=1}^{n} x\_{ij} = a\_i, \ i = 1, \ldots, m \tag{30}$$

$$\sum\_{i=1}^{m} x\_{ij} \le b\_{j\prime} \; j = 1, \dots, n \; \text{n} \tag{31}$$

$$\text{tr}\_{\text{ij}} \ge 0, \text{ } i = 1, \dots, m, \text{ } j = 1, \dots, n \tag{32}$$

#### *4.1. Container Transportation Problem*

The Container Transportation Problem is a special case of Hitchcock's transportation problem, where we assume that materials from suppliers to customers are transported only in containers of a certain capacity. Instead of rates per unit of material transported, there are prices per container transported, being fixed even if the container is not completely full.

From the previous three possibilities for the sum of all stocks and the sum of all requirement relations, the case of the stocks being sufficient to meet all the requirements is given here.

Assume that *K* is the capacity of the container and *yij* gives the number of containers needed for the quantity of material *xij*. Obviously, *yij* must be integers, the last container to reach the quantity *xij* need not be full.

Then, the container transportation problem for all requirements met can be formulated as the following model:

$$z = \sum\_{i=1}^{m} \sum\_{j=1}^{n} c\_{ij} y\_{ij} \to \min \tag{33}$$

subject to

$$\sum\_{j=1}^{n} x\_{ij} \le a\_{i\prime} \text{ i } i = 1 \text{ } \dots \text{ } \nmid m \tag{34}$$

$$\sum\_{i=1}^{m} x\_{ij} = b\_{j}, \ j = 1, \ldots, n$$

$$\text{If } \mathbf{x}\_{i\mathbf{j}} \le \mathbf{K} \\ \mathbf{y}\_{i\mathbf{j}}, \text{ } i = 1, \dots, m, \text{ } j = 1, \dots, n \text{ } \tag{36}$$

$$\{y\_{ij} \in \mathbb{Z}\_+, \ i = 1, \ldots, m, \ j = 1, \ldots, n\} \tag{37}$$

$$\mathbf{x}\_{i\bar{j}} \ge \mathbf{0}, \ i = 1, \ldots, m, \ j = 1, \ldots, n \tag{38}$$

#### *4.2. Allocation Problem*

For the transportation problem described in Section 4, it was possible to provide the required quantity by composing partial quantities from different suppliers (from different warehouses) when fulfilling the requirements.

However, in the *Allocation Problem*, it is required that the required quantity is provided from a single location so the mathematical model of the transportation problem has to be modified by [10] to account for this condition. With the same notation used for the symbols, the meaning of the decision variables *xij* is now different. They only assume binary values and *xij* = 1 if the quantity *bj* required by the *j*th customer is sourced from the *i*th supplier, if not, *xij* = 0.

If more than one customer receives the required quantity from the same supplier, the sum of their requirements must not exceed the capacity of that supplier (stock).

The model of the allocation problem with these conditions then takes the following form:

$$z = \sum\_{i=1}^{m} \sum\_{j=1}^{n} c\_{ij} \mathbf{x}\_{ij} \to \min \tag{39}$$

subject to

$$\sum\_{i=1}^{m} \mathbf{x}\_{ij} = 1, \ j = 1, \dots, n$$

$$\sum\_{j=1}^{n} b\_j \mathbf{x}\_{i\bar{j}} \le a\_{i\prime} \; i = 1, \; \dots \; \; m \tag{41}$$

$$\mathbf{x}\_{ij} \in \{0, 1\}, \ i = 1, \ldots, m, \ j = 1, \ldots, n \tag{42}$$

#### *4.3. Location Problem*

The *Location Problem* is an extension of the allocation problem [12]. For clarity, let us first summarize all the symbols used.

Consider *m* locations with capacities *ai*, *i* = 1, ... , *m* that can be used to operate warehouses supplying *n* customers with demands *bj*, *j* = 1, ... , *n*. The operation of the warehouse at the *i*th location requires a cost *fi*, *i* = 1, ... , *m* for the given period. Let *cij*, *i* = 1, ... , *m*, *j* = 1, ... , *n* be the cost of the *j*th customer being assigned to get the required quantity from the *i*th location.

The aim is to decide in which locations to operate the warehouses and to find the assignment of customers to the operated warehouses so that the value of the total cost of operating the system is minimal. Like in the allocation problem, we assume that the demands of each consumer must be covered from a single warehouse.

Therefore, the meaning of the binary decision variables *xij* is analogous to the allocation problem, *xij* = 1, if the quantity *bj* required by the *j*th customer is provided from the warehouse at the *i*th location, if not, *xij* = 0.

In addition, there are other binary decision variables *yi*, *i* = 1, ... , *m*, where *yi* = 1 means that the warehouse at the *i*th location will be operated and, if *yi* = 0, it will not be operated there.

The model of the location problem with these conditions has the following form:

$$z = \sum\_{i=1}^{m} \sum\_{j=1}^{n} c\_{ij} x\_{ij} + \sum\_{i=1}^{m} f\_i y\_i \to \min \tag{43}$$

subject to

$$\sum\_{i=1}^{m} x\_{ij} = 1, \ j = 1, \dots, n$$

$$x\_{i\bar{j}} \le y\_{i\nu} \text{ i } = 1, \dots, m, \ j = 1, \dots, n \tag{45}$$

$$\sum\_{j=1}^{n} b\_{j} \chi\_{ij} \le a\_{i\prime} \text{ i } = 1, \dots, m \text{ }\tag{46}$$

$$\mathbf{x}\_{i\bar{j}} \in \{0, 1\}, \ i = 1, \ldots, m, \ j = 1, \ldots, n \tag{47}$$

$$y\_i \in \{0, 1\}, \ i = 1, \ldots, m \tag{48}$$

As in the allocation problem, the condition (44) means that each customer takes the entire requested quantity from a single location, the condition (46) monitors the nonoverstocking of individual locations by customers receiving the requested quantity from the same location.

Let us have a look at condition (45). The left and right sides are binary variables with the inequality satisfied for the combinations 0 ≤ 0, 0 ≤ 1 and 1 ≤ 1, but not for the combination 1 ≤ 0. This ensures that no customer can get anything from a location where the warehouse will not be operated.

It is clear that, for all combinations of indices *i* and *j*, (45) represents a system of *mn* conditions. Expressing this for the values of the indices *j*, we get:

$$\begin{aligned} \mathbf{x}\_{i1} \le y\_{i\nu} \ i = 1\_{\nu} \dots i\_{\nu} m \\ \mathbf{x}\_{i2} \le y\_{i\nu} \ i = 1\_{\nu} \dots i\_{\nu} m \\ \vdots \\ \mathbf{x}\_{in} \le y\_{i\nu} \ i = 1\_{\nu} \dots i\_{\nu} m \end{aligned}$$

Summing up the previous inequalities, we get:

$$\mathbf{x}\_{i1} + \mathbf{x}\_{i2} + \dots + \mathbf{x}\_{in} \le y\_i + y\_i + \dots + y\_i \text{ i } = 1, \dots, m, \mathbf{y}$$

and hence *<sup>n</sup>*

$$\sum\_{j=1}^{n} x\_{ij} \le n y\_{i\prime} \; i = 1, \ldots, m \; \tag{49}$$

Equation (49) is equivalent to (45), but is simpler because it represents only *m* conditions rather than the *mn* conditions in the original expression (45).

#### *4.4. Capacitated Network Area Coverage*

Let us consider two finite sets *I* and *J*, where *I* is the set of service centers 1, 2, ... , *m*, and *J* is the set of customer locations 1, 2, . . . , *n*.

Further, *aij* = 1 means that customer location *j* is in a reachable distance to service center *i*, *aij* = 0 means that it does not satisfy it, and *wi* expresses the weights of service centers (since it is the minimization problem, the greater the weights are, the smaller the coefficient must be).

Similarly, *xi* = 1 means that service centre *i* is selected, while *xi* = 0 means that it is not selected.

Finally, *ci*, *i* ∈ *I*—capacity of service centre *i*, *bj*, *j* ∈ *J*—demand of customer location *j*, *yij* ∈ {0, 1}—customer from location *j* is assigned or is not assigned to service centre *i*.

In [13], we derived the following model for a capacitated network area coverage:

$$z = \sum\_{i \in I} w\_i x\_i \to \min \tag{50}$$

subject to

$$\forall j \in J: \sum\_{i \in I} a\_{ij} \mathbf{x}\_i \ge 1 \tag{51}$$

$$\forall j \in J: \sum\_{i \in I} a\_{ij} y\_{ij} = 1 \tag{52}$$

$$\forall i \in I: \ c\_i \mathbf{x}\_i \ge \sum\_{j \in J} a\_{ij} y\_{ij} b\_j \tag{53}$$

$$\forall i \in I: \sum\_{j \in J} y\_{ij} \le n\chi\_i \tag{54}$$

$$\forall i \in I: \ x\_i \in \{0, 1\} \tag{55}$$

$$(\forall i \in I)(\forall j \in J):\ y\_{ij} \in \{0, 1\}.\tag{56}$$

A necessary precondition for finding a solution is that the sum of all capacities is sufficient to cover all demands, i.e., ∑*<sup>m</sup> <sup>i</sup>*=<sup>1</sup> *ci* ≥ ∑*j*∈*<sup>J</sup> bj*, with each customer having a reachable distance to at least one center, i.e., ∀*<sup>j</sup>* ∈ *<sup>J</sup>* : <sup>∑</sup>*i*∈*<sup>I</sup> aij* > 0.

In [13], we then modified the previous model for the domain of telecommunication signals considering signal interference and its nonlinear version linearized as follows:

$$z = \left(\sum\_{i \in I} w\_i \mathbf{x}\_i\right) / \sum\_{i \in I} w\_i - \left(\sum\_{i \in I} \sum\_{j \in I} d\_{ij} l\_{i\bar{j}}\right) / \left(\sum\_{i \in I} \sum\_{j \in I} d\_{i\bar{j}}\right) \to \min \tag{57}$$

subject to

$$(\forall i \in I)(\forall j \in I):\ h\_{ij} \le x\_i \tag{58}$$

$$(\forall i \in I)(\forall j \in I):\ h\_{ij} \le x\_j \tag{59}$$

$$(\forall i \in I)(\forall j \in I):\ h\_{ij} \ge (x\_i + x\_j - 1)\tag{60}$$

$$(\forall i \in I)(\forall j \in I):\ h\_{ij} \in \{0, 1\}\tag{61}$$

$$\forall j \in J: \sum\_{i \in I} a\_{ij} \mathbf{x}\_i \ge 1 \tag{62}$$

$$\forall j \in J: \sum\_{i \in I} a\_{ij} y\_{ij} = 1 \tag{63}$$

$$\forall i \in I: \ c\_i \mathbf{x}\_i \ge \sum\_{j \in J} a\_{ij} y\_{ij} b\_j \tag{64}$$

$$\forall i \in I: \sum\_{j \in J} y\_{ij} \le m\chi \tag{65}$$

$$(\forall i \in I)(\forall j \in I)(i \neq j):\ d\_{ij} \ge (\mathbf{x}\_i + \mathbf{x}\_j - \mathbf{1})d\_{\min} \tag{66}$$

$$\forall i \in I: \ x\_i \in \{0, 1\} \tag{67}$$

$$(\forall i \in I)(\forall j \in J):\ y\_{ij} \in \{0, 1\}.\tag{68}$$

Another possible modification of the model is to meet the demand by composing parts of the capacities of several centers, but with a fragmentation not lower than a certain threshold. This new approach will be presented in detail in a separate paper.

#### *4.5. Transportation Problem with Supply from Primary Source*

Consider now a transportation network where, in addition to locations with warehouse stocks and customer requirements, there will also be a primary source, which can represent the location of the transported commodity or a global warehouse, and customers can be supplied both from local warehouses and directly from the primary source.

Assume the constraints and denotations from the location problem and two types of transportation equipment, one with a larger capacity *k*<sup>1</sup> from the primary source to local warehouses and a cost *n*<sup>1</sup> per 1 km of travel, and the other with a smaller capacity *k*<sup>2</sup> to customers and a cost *n*<sup>2</sup> per 1 km. Denoting the distance from the primary source to the *i*th local storage by *ei*, and the distance from the primary source to the *j*th customer by *gj*, we add binary decision variables *zj* to indicate whether the *j*th customer receives the desired quantity directly from the primary source (in the positive case *zj* = 1, otherwise *zj* = 0), the model with primary source and transportation technique information has the following form:

$$z = \sum\_{j=1}^{n} \frac{b\_j}{k\_2} g\_j n\_2 z\_j + \sum\_{i=1}^{m} \sum\_{j=1}^{n} \left(\frac{b\_j}{k\_1} c\_i n\_1 + \frac{b\_j}{k\_2} d\_{ij} n\_2\right) x\_{ij} + \sum\_{i=1}^{m} f\_i y\_i \to \min \tag{69}$$

subject to

$$z\_j + \sum\_{i=1}^{m} x\_{ij} = 1, \ j = 1, \dots, n \tag{70}$$

$$\sum\_{j=1}^{n} \mathbf{x}\_{ij} \le n y\_{i\prime} \; i = 1, \; \dots, m \; \; m \tag{71}$$

$$\sum\_{j=1}^{n} b\_j \boldsymbol{\chi}\_{ij} \le a\_{i\prime} \; i = 1, \; \dots \; \; m \tag{72}$$

$$\{x\_{ij} \in \{0, 1\}, \ i = 1, \dots, m, \ j = 1, \dots, n\} \tag{73}$$

$$\{y\_i \in \{0, 1\}, \ i = 1, \ldots, m\} \tag{74}$$

$$z\_j \in \{0, 1\}, \ j = 1, \ldots, n \tag{75}$$

The fractions in the objective function according to (69) indicate how many times the distance must be traveled for the customer to receive the required quantity. Since the fractions may have non-integer values, they must be rounded up to integers, the corresponding capacity may not be fully used for the last trip. The expression with the first summation in the objective function corresponds to the total cost of moving material from the primary source directly to customers, and the expression with the double summation in the objective function corresponds to the total cost of moving material from the primary source to local warehouses and from there on to the customers.

If, instead of the conditions of the location problem, the simpler conditions of the allocation problem were assumed (i.e., dropping the decision as to whether or not to use a location for storage), the previous model would be simplified, the conditions (71) and (74) would be dropped and *yi* would be omitted in the last term of the objective function (i.e., the fixed costs of all locations would be included).

#### *4.6. Crop Problem*

In plant production, an important task is to find a method of sowing the land with agricultural crops (cultures) in such a way that, given the expected yield of crops on the land and the profit from the sale of individual crops, the total profit is maximised.

Assume the following notation:

*pi*, *i* = 1, . . . , *m* ... grounds,

*ri*, *j* = 1, . . . , *m* ... area of grounds (plays the role of available capacities),

*kj*, *j* = 1, . . . , *n* ... agricultural crops (cultures),

*cij*, *i* = 1, . . . , *m*, *j* = 1, . . . , *n* ... profit from 1 ha of ground *pi*, sown with culture *kj*

*xij*, *i* = 1, . . . , *m*, *j* = 1, . . . , *n* ...number of hectares of ground *pi* sown with crop *kj*.

The mathematical model of the Crop Problem is similar to the basic version of the transportation problem with unbalanced capacities (the ground areas may not be fully used), but it lacks a set of constraints corresponding to the fulfillment of the requirements with the difference that the problem being a maximization one. It takes the following form:

$$z = \sum\_{i=1}^{m} \sum\_{j=1}^{n} c\_{ij} x\_{ij} \to \max \tag{76}$$

subject to

$$\sum\_{j=1}^{n} x\_{ij} \le r\_{i\prime} \text{ } i = 1, \dots, m \text{ } \tag{77}$$

$$\mathbf{x}\_{ij} \ge \mathbf{0}, \ i = 1, \ldots, m, \ j = 1, \ldots, n \tag{78}$$

Equation (77) expresses the use of grounds, which corresponds to the drawdown of supplier stocks in the distribution problem.

If we required each crop *j* to be sown on some minimum area *dj*, then the problem would become an example of the maximization version of the *generalized distribution problem* and the model would be modified as follows:

$$z = \sum\_{i=1}^{m} \sum\_{j=1}^{n} c\_{ij} x\_{ij} \to \max \tag{79}$$

subject to

$$\sum\_{j=1}^{n} x\_{ij} \le r\_i, \ i = 1, \dots, r, m \tag{80}$$

$$\sum\_{i=1}^{m} \mathbb{1}\_{ij} \ge d\_{j\prime} \text{ } j = 1, \dots, m \text{ }\tag{81}$$

$$\mathbf{x}\_{ij} \ge \mathbf{0}, \ i = 1, \dots, m, \ j = 1, \dots, n \tag{82}$$

#### **5. Scheduling Problems**

Scheduling problems are numerous and varied. They arise in diverse areas such as flexible manufacturing systems, production planning, computer design, logistics, timetabling, communication, etc. [14].

Here we focus on one of them, the Permutation Flow Shop Scheduling Problem (PFSSP), which, like the Assignment problem and the Traveling Salesman Problem, is permutational in nature.

It can be briefly described as follows: There are a set of *m* machines (processors) and a set of *n* jobs. Each job comprises a set of *m* operations which must be done on different machines. All jobs have the same processing operation order when passing through the machines. There are no precedence constraints among operations of different jobs. Operations cannot be interrupted, and each machine can process only one operation at a time. The problem is to find the job sequences on the machines which minimizes the makespan (i.e., the maximum completion times of all operations).

#### *Mathematical Model of PFSSP*

Consider three finite sets *J*, *M*, *O* where *J* is a set of jobs 1, ... , *n*, *M* is a set of machines 1, . . . , *m*, and *O* is a set of operations 1, . . . , *m*.

Denote:


Define the following decision variables:

∀*i*, *j* ∈ *J* : *xij* = 1, if job *j* is assigned to the *i*th position in the permutation (*Ji* = *j*) 0, otherwise (83)

> Figure 1 illustrates the use of the variables *vik* and *wik* on an example with 5 jobs and 3 machines.

**Figure 1.** Meaning of the variables *vik* and *wik*.

From Figure 1, we can draw some more general conclusions:


All verbal conclusions are expressed formally by the following system of equations:

$$\forall i \in J: \sum\_{j=1}^{n} x\_{ij} \quad = \quad 1 \tag{84}$$

$$\forall j \in J: \sum\_{i=1}^{n} x\_{ij} \quad = \quad 1 \tag{85}$$

$$\forall k \in M - \{m\}: w\_{1k} = 0 \tag{86}$$

$$\forall k \in M - \{1\}: v\_{1k} = \sum\_{r=1}^{k-1} \sum\_{i=1}^{n} p\_{ir} \mathbf{x}\_{1i} \tag{87}$$

$$\begin{aligned} \left( \forall i \in J - \{ n \} \right) \left( \forall k \in M - \{ m \} \right) &: \\ \left. v\_{i+1,k} + \sum\_{j=1}^{n} p\_{jk} \mathbf{x}\_{i+1,j} + w\_{i+1,k} \right. & \left. = \left. w\_{ik} + \sum\_{j=1}^{n} p\_{j,k+1} \mathbf{x}\_{ij} + v\_{i+1,k+1} \right. \end{aligned} \tag{88}$$

$$\mathcal{C}\_{\text{max}} \quad = \sum\_{i=1}^{n} (v\_{im} + \sum\_{j=1}^{n} p\_{jm} \mathcal{x}\_{ij}) \tag{89}$$

#### **6. Computational Results**

From the above problems, we select two, TSP and PFSSP, that are NP-hard in the decision versions [69,70].

To give an idea of the cardinality of the search space of these permutation problems, we present a few factorials as follows:

10! = 3628800 ≈ 3.6 × <sup>10</sup><sup>6</sup> 20! = 2432902008176640000 ≈ 2.4 × <sup>10</sup><sup>18</sup> 30! = 265252859812191058636308480000000 ≈ 2.6 × 1032 40! = 815915283247897734345611269596115894272000000000 ≈ 8.1 × 1047 50! ≈ <sup>3</sup> × 1064

... 100! ≈ 9.3 × 10157

The traditional approaches to such problems are based on computations using heuristic methods [71,72] for large instances such as genetic algorithms, simulated annealing, tabu search, differential evolution [73], firefly algorithm, particle swarm optimization, and ant colony optimization. Then, statistical tests are applied to examine at a certain significance level (e.g., *α* = 0.05), to what extent the mean value of the results obtained by different methods and different settings of their parameters at a larger number of runs is the same or different (and, therefore, one of the methods gives better results). For the *t*-test, we assume that the sets of values have a normal distribution. However, this assumption may be false, and then one of the non-parametric tests, such as the Wilcoxon test, must be used.

Since, given the validity of the No Free Lunch Theorem [74,75], one should not expect a general conclusion that any of the heuristics for each problem instance gives better results than other heuristics.

In this paper, we do not explore heuristics using instead a mixed integer programming model with software tools built as solvers in the GAMS environment [76,77] to find an exact solution by deterministic computation.

Statistical evaluations are, therefore, meaningless here. What can be said, however, is that the power of this software today is considerably greater than it was 20 years ago, when, in our experience, for a problem with a complexity of O(20!), the system ended up with a runtime error and the message "insufficient space to update U-factor ...". The performance of GAMS has been steadily increasing over the years, although the source code of the solvers is not freely available, from [78] it can at least be seen that it includes, among others, CPLEX, GUROBI, Lindo, and the results of the work of academic departments of Princeton University, Stanford University, and Zuse Institute Berlin. Today, GAMS calculates the exact solution for PFSSP with 20 jobs on a laptop with Intel(R) Core(TM) i5-10210U CPU @ 1.60 GHz 2.11 GHz processor, 8 GB operational memory and 64-bit operating system in less than 3 min, as shown in the following subsection.

Of course, with a computer of better technical parameters for the same time limit we get results for larger instances of the problem, but it seems to be better to use a heuristic beyond this boundary, e.g., our GA 'war elimination' modification [79].

Since PFSSP has a more complex model than TSP, we start with it and include its complete GAMS code.

#### *6.1. PFSSP Computational Results*

For PFSSP with 10 jobs, 6 machines, processing times from the TABLE section (it corresponds to the first benchmark in Table 1) and the model given by Equations (84)–(89), the program code in GAMS can be, e.g., as follows:

```
* Permutation flow shop scheduling problem
$TITLE permutation flow shop scheduling problem
$OFFSYMXREF
$OFFUELLIST
$OFFUELXREF
```

```
OPTION ITERLIM=200,000
* ITERLIM number of iterations
OPTION OPTCR=0.00001
*OPTION OPTCR=0.001
* OPTCR stopping in MIP in case the best solution is within the limits
* 100*OPTCR% of the global~extreme
* section defining indexes
SETS
  I jobs /1*10/
```

```
K machines /1*6/;
```

```
ALIAS(I,J);
* J - position of the job in the permutation
  ALIAS(K,R);
* input data section
PARAMETERS
  N,M;
  N=CARD(I);
  M=CARD(K);
TABLE P(I,K)
      123456
  1 333 991 996 123 145 234
  2 333 111 663 456 785 532
  3 252 222 222 789 214 586
  4 222 204 114 876 752 532
  5 255 477 123 543 143 142
  6 555 566 456 210 698 573
  7 558 899 789 124 532 12
  8 888 965 876 537 145 14
  9 889 588 543 854 247 527
 10 999 889 210 632 451 856;
*variables section (decision variables and objective function)
VARIABLES
  X(I,J) is 1 if job j is assigned to position i in the permutation, 0, otherwise
  V(I,K) waiting time on machine k before the start of job i in the permutation
  W(I,K) waiting time of job i in the permutation after finishing processing
* on machine k, while machine k+1 becomes free
  Cmax total processing time for all tasks (makespan);
BINARY VARIABLE X;
NONNEGATIVE VARIABLE V;
NONNEGATIVE VARIABLE W;
*section describing the system of (in)equalities
EQUATIONS
  EQ1(I)
  EQ2(J)
  EQ3(K)
  EQ4(K)
  EQ5(I,K)
  OBJFCE(K);
  EQ1(I) .. SUM(J,X(I,J)) =E= 1;
  EQ2(J) .. SUM(I,X(I,J)) =E= 1;
  EQ3(K)$(ORD(K) LE (M-1)) .. W(''1'',K) =E= 0;
  EQ4(K)$(ORD(K) GE 2)
      .. V(''1'',K) =E= SUM((R,I)$(ORD(R) LE (ORD(K)-1)),P(I,R)*X(''1'',I));
  EQ5(I,K)$((ORD(I) LE (N-1)) AND (ORD(K) LE (M-1)))
    .. V(I+1,K)+SUM(J,P(J,K)*X(I+1,J))+W(I+1,K) =E=
      W(I,K)+SUM(J,P(J,K+1)*X(I,J))+V(I+1,K+1);
  OBJFCE(K)$(ORD(K) EQ M) .. Cmax =E= SUM(I,V(I,K)+SUM(J,P(J,K)*X(I,J)));
*description of the model, running the solver, and displaying the results
MODEL FLOWSHOP /ALL/;
SOLVE FLOWSHOP USING MIP MINIMIZING Cmax;
DISPLAY X.L, V.L, W.L, Cmax.L;
```
The ability to compute optimal solutions was checked using standard benchmarks from OR-Library (OR = Operations Research) accessible at Brunel University London [80], originally described in [81]. The computational results are summarised in Table 1. For

instances with 30 or more jobs, GAMS does not find the optimal solution in the 1000 s time limit, but only a "close" approximation, which, however, differs by less than 10% even for the last instance 75 × 20, where the optimal value is unknown due to the huge size of the search space and is only estimated by the interval. For such cases, we at least suggest a solution method using the genetic algorithm [82].

**Table 1.** GAMS computational results (10 × 6 corresponds to 10 jobs and 6 machines, etc.; t-l-e = time limit exceeded).


To do this, we will need a model that builds an appropriate schedule for the permutation. The genetic algorithm will then select a promising part of the search space of permutations in which a good approximation of the optimum can be found in a reasonable amount of time.

If we have processing times *pij* for job *i* on machine *j*, and a job permutation *J*1, *J*2, ... , *Jn*, then we can calculate the completion times *CJi*,*<sup>j</sup>* as follows:

$$C\_{l\_1,1} = \\_p\_{l\_1,1} \tag{90}$$

$$\forall i \in I - \{1\}: \mathbb{C}\_{I\_i, 1} \quad = \quad \mathbb{C}\_{I\_{i-1}, 1} + p\_{I\_i, 1} \tag{91}$$

$$\forall k \in M - \{1\}: \mathbb{C}\_{I\_1, k} = \begin{array}{c} \mathbb{C}\_{I\_1, k-1} + p\_{I\_1, k} \end{array} \tag{92}$$

$$(\forall i \in I - \{1\})(\forall k \in M - \{1\}): \mathbb{C}\_{l\_i k} = \max\left\{\mathbb{C}\_{l\_{i-1} k}, \mathbb{C}\_{l\_i k - 1}\right\} + p\_{l\_i k} \tag{93}$$

$$\mathbb{C}\_{\text{max}} \quad = \quad \mathbb{C}\_{\text{l}\_{\text{u}}, \text{m}} \tag{94}$$

As the genetic algorithms [79] are well known, we only summarise parameter settings and describe only the problem-specific operators in more detail.

The fitness function is inversely proportional to the makespan, the smaller the makespan, the higher the value of the fitness function.

The number of individuals in the population was set to 50 and the number of iterations to 10*n*2. The initial population was generated randomly, and the parents for the crossover operation were determined by binary tournament selection.

As to the *crossover* operation, we cannot use the traditional two-point crossover, because it would lead to infeasible solutions. If we change the middle parts of the parent chromosomes *P*<sup>1</sup> and *P*<sup>2</sup> in Figure 2, we will obtain offspring (10,5,2,6,10,4,1,6,3,1) and (5,8,4,7,8,9,2,6,3,1) that correspond to no permutations, because some jobs are duplicated or omitted. We used what is called *crossover in a partially mapped representation* where the genes in the middle part of one chromosome are ordered in its offspring by their occurrence in the second parent chromosome.


**Figure 2.** Modified two-point crossover.

From the possible mutation operations, we have selected the *shift mutation*, which removes a value at one position and puts it at another position), see Figure 3.

**Figure 3.** Shift mutation.

The offspring of the parents replaced two randomly selected individuals with belowaverage fitness function values.

With the above parameter settings, the results for the last 4 instances in Table 1 were obtained as shown in Table 2. The average values from 30 runs are presented here, as well as the best values obtained from them, which for these large instances are better than the values obtained from GAMS when the 1000 s timeout expires. All these results were achieved in less than 10 s because of the small number of iterations of the genetic algorithm.

**Table 2.** GA computational results—average and the best result from 30 runs, optimum.


#### *6.2. TSP Implementation in GAMS*

In describing the source code in GAMS and verifying its computational abilities, we use three benchmarks from the TSPlib library [83] with 24, 52, and 100 cities, or positions in the map given by coordinates.

The following code is written for the gr24.tsp benchmark. Since the adjacency matrix is symmetric, only the data of the lower triangular matrix are entered with the remaining data calculated. The EQUATIONS section is a rewrite of the TSP model and its Equations (5)–(8). The *xij* binary domain, corresponding to Equation (9), is given by the declaration that precedes this section.

```
$TITLE Travelling Salesman Problem
OPTION ITERLIM=10000000;
OPTION OPTCR=0;
```

```
SETS
  I /1*24/;
```

```
ALIAS (I,J);
PARAMETERS
 N;
 N=CARD(I);
TABLE C(I,J) adjacency matrix
     1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
1 0
2 257 0
3 87 196 0
4 91 228 158 0
5 150 112 96 120 0
6 80 196 88 77 63 0
7 130 167 59 101 56 25 0
8 134 154 63 105 34 29 22 0
9 243 209 286 159 190 216 229 225 0
10 185 86 124 156 40 124 95 82 207 0
11 214 223 49 185 123 115 86 90 313 151 0
12 70 191 121 27 83 47 64 68 173 119 148 0
13 272 180 315 188 193 245 258 228 29 159 342 209 0
14 219 83 172 149 79 139 134 112 126 62 199 153 97 0
15 293 50 232 264 148 232 203 190 248 122 259 227 219 134 0
16 54 219 92 82 119 31 43 58 238 147 84 53 267 170 255 0
17 211 74 81 182 105 150 121 108 310 37 160 145 196 99 125 173 0
18 290 139 98 261 144 176 164 136 389 116 147 224 275 178 154 190 79 0
19 268 53 138 239 123 207 178 165 367 86 187 202 227 130 68 230 57 86 0
20 261 43 200 232 98 200 171 131 166 90 227 195 137 69 82 223 90 176 90 0
21 175 128 76 146 32 76 47 30 222 56 103 109 225 104 164 99 57 112 114 134 0
22 250 99 89 221 105 189 160 147 349 76 138 184 235 138 114 212 39 40 46 136 96 0
23 192 228 235 108 119 165 178 154 71 136 262 110 74 96 264 187 182 261 239 165 151 221 0
24 121 142 99 84 35 29 42 36 220 70 126 55 249 104 178 60 96 175 153 146 47 135 169 0;
SET C2(I,J);
C2(I,J)$(NOT SAMEAS(I,J)) = yes;
C(C2(I,J)) = MAX(C(I,J),C(J,I));
VARIABLES
 X(I,J)
 delta(I)
 Z;
BINARY VARIABLE X(I,J);
EQUATIONS
 EQ1(J) each city is entered exactly once
 EQ2(I) each city is left exactly once
 EQ3(I,J) subtour elimination constraints
 EQ4 objective function;
 EQ1(J) .. SUM(I,X(I,J)$(ORD(I) NE ORD(J))) =E= 1;
 EQ2(I) .. SUM(J,X(I,J)$(ORD(I) NE ORD(J))) =E= 1;
 EQ3(I,J)$((ORD(I) GE 2) AND (ORD(J) GE 2) AND (ORD(I) NE ORD(J)))
        .. delta(I)-delta(J)+N*X(I,J) =L= N-1;
 EQ4 .. Z =E= SUM((I,J),C(I,J)*X(I,J));
MODEL TSP/ALL/;
SOLVE TSP USING MIP MINIMIZING Z;
DISPLAY X.L, Z.L;
```
The total length of the route is 1272, the decision variables *xi*,*<sup>j</sup>* have a value of 1 in the following order: *x*1,16, *x*16,11, *x*11,3, *x*3,7, *x*7,6, *x*6,24, *x*24,8, *x*8,21, *x*21,5, *x*5,10, *x*10,17, *x*17,22, *x*22,18, *x*18,19, *x*19,15, *x*15,2, *x*2,20, *x*20,14, *x*14,13, *x*13,9, *x*9,23, *x*23,4, *x*4,12, *x*12,1, and, thus, the circuitous route passes through cities 1, 16, 11, 3, 7, 6, 24, 8, 21, 5, 10, 17, 22, 18, 19, 15, 2, 20, 14, 13, 9, 23, 4, 12, 1, which is in agreement with the published result for the gr24 benchmark. The calculation time was 0.36 s.

The data for the berlin52.tsp benchmark are entered differently, namely as a matrix with 52 rows and 2 columns, where the 1st column is the value of the *x*-coordinate and the 2nd column is the value of the *y*-coordinate. The adjacency matrix in this case is obtained by calculating the Euclidean distances between all pairs of positions. This part of the code takes the following form, the rest is the same as in the previous code.

SETS I /1\*52/; ALIAS (I,J); PARAMETERS N;

N=CARD(I);


PARAMETERS C(I,J);

SET C2(I,J); C2(I,J)\$(NOT SAMEAS(I,J)) = yes; C(C2(I,J)) = ROUND(SQRT(SQR(XY(I,'1')-XY(J,'1'))+SQR(XY(I,'2')-XY(J,'2'))));

The total length of the route is 7542, the calculation time was 1.45 s and the circuitous route passes through positions 1, 49, 32, 45, 19, 41, 8, 9, 10, 43, 33, 51, 11, 52, 14, 13, 47, 26, 27, 28, 12, 25, 4, 6, 15, 5, 24, 48, 38, 37, 40, 39, 36, 35, 34, 44, 46, 16, 29, 50, 20, 23, 30, 2, 7, 42, 21, 17, 3, 18, 31, 22, 1. The optimal route can be seen in Figure 4.

**Figure 4.** The optimal route for the berlin52.tsp benchmark.

Finally, GAMS for the kroA100.tsp benchmark stopped the computation at 1000.02 s by exceeding the time limit, but the intermediate result of the path length 21282 and its traversal through positions 1, 47, 93, 28, 67, 58, 61, 51, 87, 25, 81, 69, 64, 40, 54, 2, 44, 50, 73, 68, 85, 82, 95, 13, 76, 33, 37, 5, 52, 78, 96, 39, 30, 48, 100, 41, 71, 14, 3, 43, 46, 29, 34, 83, 55, 7, 9, 57, 20, 12, 27, 86, 35, 62, 60, 77, 23, 98, 91, 45, 32, 11, 15, 17, 59, 74, 21, 72, 10, 84, 36, 99, 38, 24, 18, 79, 53, 88, 16, 94, 22, 70, 66, 26, 65, 4, 97, 56, 80, 31, 89, 42, 8, 92, 75, 19, 90, 49, 6, 63, 1 corresponds to the known optimal solution for this benchmark. The optimal route is in Figure 5.

**Figure 5.** The optimal route for the kroA100.tsp benchmark.

For instances with more than 100 positions, it would be necessary to search for an approximation of the optimum using one of the heuristic methods.

One of the first was the use of the so-called *Lin-2-Opt change operator* [8], see Figure 6. Here, two elements are added to the permutation of *n* cities to visit (into positions 0 and *n* + 1), and then the starting city is assigned to those positions to simulate a cyclic tour. Two 'edges' (pairs of neighbouring elements in permutation) are randomly chosen ((*p*1, *p*2) and

(*q*1, *q*2) say), the inner elements *p*2, *q*<sup>1</sup> are swapped and the elements between *p*<sup>2</sup> and *q*<sup>1</sup> are reversed.

**Figure 6.** Lin-2-Opt change neighbourhood operation,

Positions 0 and *n* + 1 with the fixed value of the starting city can also be used to expand the individuals in the population of the genetic algorithm and then apply the operations presented for the PFSSP.

However, we no longer investigate this for extremely large instances of the TSP problem because heuristics do not guarantee finding an optimal solution, which often is not even known here, and the aim was to find bounds for which we still obtain the precise solution in reasonable time using a 'normal' computer. In the case of GAMS, this bound is an instance with 100 cities.

#### *6.3. Data, Changes in Time, Uncertainty*

Data from OR-Library and TSPLIB are related to a specific point in time, in reality they may change over time or may not be completely known.

A more general case of the Travelling Salesman Problem is the *Canadian Traveller Problem* (CTP) [84,85]. Here, the distance matrix may change over time due to the occurrence of events that make some parts of the route inaccessible so that an adaptive strategy must be found. These events are random in nature, which corresponds to the problems of robot navigation in environments where the distribution of obstacles is only discovered as the robot moves through the environment; moreover, the obstacles may move, and thus the locations of potential collisions change dynamically.

In transport tasks, the values of some parameters can change over time, the fuel price is not constant, and the vehicle consumption can only be estimated because it can change according to the traffic situation and the season, which will affect, e.g., the calculation of the objective function (69). Similarly, in the crop problem, we can only estimate crop yields.

In location-based tasks, the problem may arise of adding another center to an existing network of centers to improve the coverage of an area. An example might be an expansion of the existing supermarket network of a chain store. Here it is suggested to use one of the properties of the Voronoi diagram [86,87], a data structure known from computational geometry: Assume a Voronoi diagram with its sites represented by the current centers. The point *q* is the vertex of the Voronoi diagram if and only if the *largest empty circle C*(*q*) contains three (or more in a degenerate case) sites on its boundary and none inside. Among these circles, we determine the one with the largest diameter, and its center is then the optimal position for the location of the new center, see Figure 7.

In fact, the calculated position may not be available, the cost of building here may be too high, thus a suitable nearby location must be found, or the center of one of the other empty circles must be chosen in descending order of diameters.

**Figure 7.** Finding a new location using the largest empty circles.

Another issue is the amount of inventory in logistics operations. The stock changes over time according to demand and needs to be replenished accordingly. We speak about *inventory management* and the resulting sustainability [88–91]. However, demands are stochastic in nature and, in addition, inventory management must take into account the cost of maintaining inventory and losses from premature depletion of inventory for undelivered goods.

In [92], a new mathematical model is derived, the properties of the profit function are proved, and the profitability in a two-channel production system considering carbon emissions and green technology is numerically verified on specific data.

While artificial neural networks (ANN) have very little application in combinatorial optimization, their main use is in cluster analysis, pattern recognition, image processing and prediction, in [93] the authors present an efficient implementation of ANNs in an inventory management model under uncertainty and inflation.

In [94] the unreliability of the supply chain and methods to eliminate this unreliability are explored, and the required mathematical equations are derived and verified by numerical experiments, including sensitivity analysis.

However, all these aspects are beyond the scope of this paper and can be the subject of separate texts as also evidenced by the papers mentioned.

#### **7. Conclusions**

This paper studies the assignment problem and its modifications with logistics applications, in routing, distribution, and scheduling tasks. Its first contribution is the correlation of the problem models, which are often distant in nature and time complexity.

It has also shown how the described models can be directly transferred to the GAMS environment. NP-hard Permutation Flow Shop Scheduling Problem (PFSSP) and the Travelling Salesman Problem are used to show that the optimal solution can be determined in the available time of a few minutes for instances with 20 jobs on 10 machines in the case of PFSSP, and for 100 cities in the case of TSP.

Previously, these boundaries were inaccessible with mixed integer programming solvers, but with the new version of GAMS, they have been significantly extended. This of course means first to build the appropriate model (and this is not always a simple matter, as the informal derivation of the PFSSP model in Section 5 showed) and then, for instance, for benchmark libraries (e.g., OR-Library or TSPLIB), to search individually for the appropriate bounds. The findings from PFSSP and TSP are not isolated examples of the successful application of GAMS in solving large instances of optimization versions of NPcomplete/NP-hard problems. We have already validated it in [13] in solving the covering problem with matrices of hundreds of thousands of elements, and more recently in solving the Steiner problem in graphs in [95], where first using the terminology of network flows

a mixed integer programming model was derived, then modified for GAMS, and finally exact results for a representative class of benchmarks from OR-Library were obtained.

Another goal of this paper was to introduce code generation in GAMS on non-trivial tasks because in the manuals [76,77] we can find only a description of individual elements of this tool, but not the codes of complete task models. In MATLAB, running the computation of an optimization program means writing just a single command (intlinprog or linprog with the appropriate parameters). Similarly, when solving differential equations, e.g., to calculate the differential equation *y* (*x*) = 4*xy* + *x*<sup>3</sup> with initial condition *y*(4) = 2, it is enough to enter dsolve('Dy=4\*x\*y+x∧3','y(4)=2','x'). In MATHEMATICA, too, to obtain the impulse function of the system described by a differential equation, it is enough to rewrite it in the form of a Laplace transfer and use a single command InverseLaplaceTransform. In contrast, the code notation in GAMS is similar to code in programming languages with the definition of constants, the declaration of variables, and the body of the program. Again, there are assignment statements, conditional statements, and loop statements. For example, the binary values of the reachability matrix **A** from the distance matrix **D** and the defined reachable distance threshold *Dmax* are determined in GAMS as follows:

LOOP(I,

```
LOOP(J,
           IF (D(I,J) <= Dmax,
                 A(I,J)=1;
                ELSE
                 A(I,J)=0;
              );
       );
);
```
The only disadvantage of GAMS is that it has no graphical tools, and the results of the calculations are only in text form. This requires exporting them to a suitable program and postprocessing. In [13] we used MATLAB, here Figures 4 and 5 are generated in the MATHEMATICA environment.

Only where for extremely large instances of problems of exponential complexity we cannot obtain an exact solution using GAMS in a reasonable amount of time (e.g., no more than in tens of minutes), do we use one of the many heuristic methods. Given the No Free Lunch Theorem [74,75], none of them can be recommended as the best in the general case, since finding the optimal solution is not guaranteed and the result is always an approximation of the optimum, so our modification of the genetic algorithm, implemented in Java and described in more detail in [79], can be used without loss of generality.

One-point heuristics (hill climbing, tabu search, simulated annealing) in solving problems where in each iteration the neighborhood operation often generates tens of infeasible solutions and it is necessary to use a repair operator for them (here it concerns the coverage problem), and slow down the computation considerably, so in these cases we prefer, e.g., a genetic algorithm that generates only two new solutions in each iteration.

The model in Section 4.4 is original with another possible modification proposed at the end. Its model has already been built and verified on smaller-scale instances so far and will be investigated in more complex cases.

In future research, we expect to focus on the Quadratic Assignment Problem, the Vehicle Routing Problem and its solvability using GAMS, and applications in agriculture with consideration of data uncertainty using probabilistic models or fuzzy modeling, since yields can only be estimated. Although the Quadratic Assignment Problem has a non-linear objective function with quadratic terms, it can be converted to a mixed integer programming problem using Lawler's linearization [7] and the MIP solver of GAMS can be used again.

**Funding:** This research received no external funding.

**Data Availability Statement:** Not applicable.

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

#### **Abbreviations**

The following abbreviations are used in this manuscript:


#### **References**

