Optimization Method to Address Psychosocial Risks through Adaptation of the Multidimensional Knapsack Problem

: This paper presents a methodological scheme to obtain the maximum benefit in occupational health by attending to psychosocial risk factors in a company. This scheme is based on selecting an optimal subset of psychosocial risk factors, considering the departments’ budget in a company as problem constraints. This methodology can be summarized in three steps: First, psychosocial risk factors in the company are identified and weighted, applying several instruments recom-mended by business regulations. Next, a mathematical model is built using the identified psychosocial risk factors information and the company budget for risk factors attention. This model represents the psychosocial risk optimization problem as a Multidimensional Knapsack Problem (MKP). Finally, since Multidimensional Knapsack Problem is NP-hard, one simulated annealing algorithm is applied to find a near-optimal subset of factors maximizing the psychosocial risk care level. This subset is according to the budgets assigned for each of the company’s departments. The proposed methodology is detailed using a case of study, and thirty instances of the Multidimensional Knapsack Problem are tested, and the results are interpreted under psychosocial risk problems to evaluate the simulated annealing algorithm’s performance (efficiency and efficacy) in solving these optimization problems. This evaluation shows that the proposed methodology can be used for the attention of psychosocial risk factors in real companies’ cases.


Introduction
In any company, adequate attention to psychosocial risk (PSR) factors and the generation of efficient institutional policies regulating interpersonal relationships allow anticipating risks associated with a low productivity level. Looking to the future, algorithmic management, rapid changes in the structure of work and the workforce [1,2] pose an unknown situation, which possibly increased the risks of work stress. Problems such as high staff turnover, difficulties in communications, poor leadership, psycho-logical stressors, coping mechanisms, family work relationship, among other factors, constitute a barrier to the progress and evolution of organizations [3,4]. The benefits of attending to these psychosocial aspects include boosting productivity, wellbeing of their employees, moving towards a highly competitive company and improving market positioning [5,6].
In addition, in Industry 4.0, current digital tools (a new generation of sensors) and recent developments in information technologies (Big Data, Machine Learning, Artificial Intelligence, Internet-of-Things) play an essential role in the development of more competitive and globalizing companies. Furthermore, these new developments impact the workforce by obtaining advantages to improve worker motivation, implementing policies for organizational and human aspects and promoting health in the workplace [7][8][9].
A large number of quality publications in the existing literature promote workers' health. For example, the international ISO-45003 standard proposes implementing "good practices" to avoid PSRs and having a healthy work environment [10,11]. Likewise, the European Union Information Agency for Occupational Safety and Health (EU-OSHA) and the Spain National Institute for Safety and Health at Work (INSST) promoting the improvement of occupational safety and health conditions to decrease occupational hazards, work accidents, and occupational diseases [12].
Furthermore, Cox et al. [13] identify ten potentially dangerous work psychosocial characteristics, classifying them according to their relationship with the work context. In Kompier and Levi [14], the checklists of the European Foundation for the Improvement of Living and Working Conditions (Eurofound) are discussed. In [15], the Guide to the Actions of the Labor and Social Security Inspectorate on Psychosocial Risks points out the advantage of quickly, validly, and reliably information collecting using questionnaires complying with established requirements such as reliability and validity. In Sebastián et al. [16], several treatment recommendations and ergonomic and psychosocial evaluations are given.
The attention of occupational PSR factors has led to the developing and compliance of regulations that consider that an unfavorable work environment can cause physical and mental illnesses [10]. These regulations aim to improve work environments so that workers' activities develop favorably.
Concerning the instruments identifying PSR factors, in Spain, both the Psychosocial Evaluation Procedure called Psychosocial Factors (FPSICO, for its acronym in Spanish) of the National Institute of Occupational Safety and Health (INSHT, for its acronym in Spanish) and the Copenhagen Psychosocial Questionnaire (CoPsoQ) of the Union Institute of Labor, Environment, and Health (ISTAS) contain questionnaires for diverse types of companies [17,18]. Furthermore, the European Working Conditions Surveys (EWCS) [19] develop instruments to identify the employment situation, such as the Labor Day planning and duration, the work organization, training and learning, and physical and psychosocial risk factors. On the other hand, the Official Mexican Standard NOM-035 [10] establishes the elements to identify, analyze and prevent PSR factors, and to promote a favorable organizational environment in the workplace.
In particular, in Leka et al. [20], the policy context for managing work-related PSRs in the European Union (EU) is discussed. They highlight the importance of properly managing these risks, with the financial benefits for workers' insurance premium payments. In Sureda [21], PSR factors in one public hospital are evaluated with the FPSICO procedure, composed of 75 ten-level Likert items and seven factors: mental workload, temporal autonomy, job content, supervision-participation, role definition, worker's interest, and personal relationships.
In light of these works, it can be seen that PSR assessment is becoming increasingly important for research on occupational safety and health. However, in [22], it is pointed out that even though several sources guide the identification of PSRs, the probability and statistical evaluation of PSRs remain poorly studied. For example, in [23], a significant correlation between work resources and work-related stress symptoms is shown.
As in other disciplines of knowledge, several propositions incorporate computational tools and artificial intelligence techniques to improve risk attention and occupational health. In Khakzad [24], a heuristic is applied for using imprecise probabilities in a Bayesian network for system safety assessment under uncertainty. In Han et al. [25], the Random Forest classification heuristic is implemented to select the optimal feature set workrelated stress measured by physiological signals of electrocardiograms (ECG) and respiration (RSP). The detection of work stress was carried out with the support of a wearable device, concluding that excessive stress decreases work efficiency and leads to negative emotions and various diseases.
In one particular case, job rotation is a successful strategy used in manufacturing systems. Job rotation helps prevent musculoskeletal disorders, eliminates boredom, and increases job satisfaction and morale, providing benefits to both workers and the company. Asensio-Cuesta et al. [26] developed a genetic algorithm for this problem, and Song et al. [27] implemented a heuristic using factors related to workload ergonomic evaluation data generated in the workplace. On the other hand, Yan et al. [28] developed data processing algorithms with real-time warning indicators for evaluating and signaling risk positions through a smartphone.
Concerning optimization processes, the Multi-dimensional Knapsack Problem (MKP) has been adapted to diverse problems requiring the selection of elements with capacity restrictions in containers. MKP has been applied to solve various practical problems in the industry, such as resource allocation, transportation [29], and production planning [30,31]. In [32], a survey of the multiple applications of 0-1 MKP is presented, emphasizing that the most popular approaches are based on metaheuristics, highlighting population-based strategies such as Genetic Algorithm (GA), Particle Swarm Optimization (PSO), and Ant Colony Optimization (ACO). Most of these procedures are evaluated using various benchmarks to demonstrate their effectiveness. Furthermore, Drake et al. [33] propose a hyper-heuristic to select the arguments of crossover operators to solve a group of instances obtained from three well-known benchmark libraries (OR-Lib, SAC-94, and that presented by Glover and Kochenberger). In [34], an approach was proposed in which a greedy algorithm is first used to reduce the search space, and then a linear programming algorithm is used to obtain the optimal values of some instances of the OR-Lib. In Dzalbs et al. [35], a variant for population selection in a GA using several instances of two benchmark libraries (OR-Lib and Glover and Kochenberg) is presented. Finally, in [36,37], both ACO and PSO are used to solve several SAC-94 instances, achieving the optimal value in all of them.
From the existing occupational risk attention literature, and to the best of our knowledge: (1) There is no evidence that some methodology uses an optimization technique generating a subset of risk factors to be included in a plan for PSR factors' attention.
(2) There are no proposals considering a budget assigned to each area or department of the company for making decisions about the prioritization of occupational risk care.
This article proposes a new methodology to find the combination of PSR factors that maximizes the company's attention level, considering the company's budgets for each of its departments. The risk factors are first identified using the questionnaires recommended by the applicable regulations. Then, the problem of optimizing PSR is modeled as a Multiple Knapsack Problem (MKP). Finally, an optimal solution for this problem is found using an adaptation of the simulated annealing (SA) algorithm.

Materials and Methods
This section describes the materials and methods used to carry out the mapping scheme to build an MKP from a PSR optimization problem, and the adaptation of the SA algorithm to find a solution.

Psychosocial Risk Factors
According to the 1984 Joint International Labour Office (ILO)/World Health Organization (WHO) Committee on Occupational Health, the PSR are "the interactions between and among work environment, job content, organizational conditions, and workers' capacities, needs, culture, personal extra-job considerations that may, through perceptions and experience, influence health, work performance, and job satisfaction" [38]. 2014/2015 European Healthy Work Campaign asserts that "psychosocial factors are linked not only to health outcomes but also to performance-related outcomes such as absenteeism, workability and especially job satisfaction" [39]. For the determination of PSR factors, multiple models and variables have been proposed, and a classification of them is needed. The related literature shows the dimensions coinciding in most of the models [4]:

•
Psychosocial stressors: These imply stressful worker demands such as the workload, autonomy offered by the position, support, or harassment from the boss.

•
Coping strategies: These are related to the internal processes of the individual, such as their perception and the way they react to events, among others.

•
Family work relationships: These are the extra-work aspects that affect the work environment. This relationship is a dimension little considered in the workplace.
In particular, the most common PSR factors are workload, dealing with other people, lack of resources (equipment or materials to carry out the work), lack of support (implying the perception of insufficient support, cooperation, and recognition), poor working conditions (excessive noise and heat), and other work aspects that are annoying for the staff.
Among the benefits of addressing these PSR factors are preventing and promoting a safe work environment. Additionally, more direct benefits are obtained, such as:

•
Improve the ability to respond to regulatory compliance issues.

•
Avoid incidents and thus reduce costs to companies.

•
Reduce downtime and downtime costs.

•
Reduce the cost of insurance premiums.

•
Reduce absences and employee turnover rates.

•
Recognition for having achieved an international goal (which in turn can influence clients who are concerned about their social responsibilities).
The effects of poor PSR management translate into high rates of temporary work disability for accidents and insurance policy claims. Therefore, the attention given to PSRs does not differ from the generic risk management process in any company. Both the International Standard ISO-45003 and the Mexican NOM-035 contemplate the following steps: identification of risk, evaluation of its impact/effect on the organization, determination of the best actions to follow, and finally, their application to generate maximum benefit in the organizational environment. Additionally, several countries' regulations require the introduction of procedures promoting improvements in the safety and health of workers at their workplace. In the European Union, the European Framework Directive 89/391/EEC establishes a legal framework indicating the need for a continuous risk assessment and its reduction. This implies a continuous improvement process that must be repeated within a period established in the organizational context [20,40,41].

Multidimensional Knapsack Problem
MKP is a 0-1 combinatorial optimization problem that has been extensively studied [42,43]. MKP is classified as an NP-hard problem [44]. The MKP is represented by the optimization model depicted in Equations (1) to (3) [44]. In Martello and Toth's books [45,46], there is an extensive explanation of the properties and algorithms for MKP resolution and its variants.

=
(1) Subject to Є 0,1 j = 1, …, n MKP is described as a set of n elements with gain pj > 0, j = {1, ..., n}, and m resources (backpacks) with capacities ci > 0 i = {1, ..., m}. Each j-th element consumes an amount aij > 0 of each i-th resource. Equation (1) represents the objective function, that is, to choose a subset from the n elements maximizing the total gain value. This selection must not exceed the maximum capacity of each backpack m, expressed in the set of constraints in inequalities (2). Finally, Equation (3) indicates that x variables take only binary values, i.e., xj = 1, if the j-th element is selected, and xj = 0 otherwise. An MKP graphical representation is shown in Figure 1, where there are n boxes with a profit pj represented by a specific color in Figure 1a. Each j-th box contains different shapes with a weight aij, shown in Figure 1b. When selecting the j-th box, each shape is stored in the backpack corresponding to it. Each backpack has a specific maximum weight capacity, which cannot be exceeded.

The Simulated Annealing Algorithm
Simulated Annealing (SA) is a local search algorithm. It is named simulated annealing by analogy to the physical process of annealing solids, where a solid crystalline structure is first heated and then slowly cooled until it reaches a stable state (a state of minimum energy) that is free from defects of crystallization. The SA algorithm is inspired by this thermodynamic behavior, which is compared to the search for an optimum value for discrete optimization problems [47][48][49]. Due to its simple implementation, convergence properties, and local optimal escape ability, SA has become a popular heuristic for solving complex problems in recent decades. Since many engineering, planning, and manufacturing problems can be modeled as minimizing or maximizing a cost function [50], they can be solved using SA-based approaches.
When an optimization problem is solved, in each SA iteration, the objective function improvement is valued by comparing two solutions, the current one and one randomly selected from its neighborhood. The best solutions are always accepted, while only a fraction of "not good" solutions are accepted to escape from one local optimum and continuing the search for better solutions. The probability of accepting "not good" solutions (i.e., moving to solutions that generate a worse value in the objective function) depends on the acceptance criterion of SA based on the Boltzmann distribution. By gradually decreasing the temperature to zero, the worst movements will be accepted less frequently, and the distribution of the solution associated with the inhomogeneous Markov chain modeling the algorithm behavior converges to a distribution, in which the probability is concentrated in a more significant part in the set of optimal solutions. It is necessary to specify the set of algorithm parameters [51][52][53] presented in Table 1.

Parameter
Description w0 Initial candidate solution T External cycle control parameter T0 Initial value of the control parameter α Control Final value of the control parameter Figure 2 shows the methodological stages for optimizing the attention of PSR factors. These stages are described in the following paragraphs:

Propoused Methodology
• Identify each j-th PSR factor vj and its level of required attention pj, by applying NOM035, ISTAS-CoPsoQ, or F-Psico-based questionnaire [10,16] or using any instrument for PSR prevention of company workers, and qualify them according to the scales integrated into the instrument.

•
Construct the cost matrices of the areas or departments that address the PSR factors (a minimum of four is suggested: training, communication, industrial safety with non-conformity mechanisms, and human resources with social support actions). Assign the maximum budget for each department.

•
Map the PSR optimization problem as an MKP using the previously obtained data.

•
Tune and apply the SA algorithm to solve the PSR factors optimization problem.

•
Generate reports to comply with regulations and provide the maximum benefit to the company's labor wellbeing. The first four stages are detailed in the next sub-sections. Stage five of this methodology is developed in the next section.

Identification of PSR Factors and Generation of Cost Matrices
Reference guide number III, included in the Mexican NOM-035 standard, identifies 25 dimensions, ten domains, and five categories of psychosocial risks in the workplace [10]. Additionally, NOM-035 contains instructions to determine the risk level. For this investigation, a study case in a company is shown in Figure 3, named CASE-1. This figure shows that the highest of the five risk levels occurs in work time management. Table 2 describes the four company departments addressing these risk factors and the actions that would be feasible to apply by each one of them.  Design brochures or posters that graphically remind workers to organize their time so as not to be overloaded by work Implement a nonconformity mechanism. Through a module on the company's website where the worker can comment or suggest alternatives for workload issues associated with their position or organizational climate Program actions that promote social supportpublicize strategies or support schemes both by the company and among colleagues, an example would be to facilitate a mechanism to share excess work or cover in special events The company's departments must apply budget resources, either monetary or not, to obtaining a budget consumption matrix. For CASE-1, the matrix has five rows and four columns, as shown in Table 3. This matrix contains the monetary costs incurred to attend to each of the company's risk aspects. The company's departments must also indicate the maximum budget that will be applied to address these risks. By maximizing attention to the psychosocial risk-weighted elements subject to departments' budget, decision makers are provided with essential information for organizational wellbeing. The variables used to build the PSR optimization model are obtained from the first two stages of the proposed methodology. Each vj, j = {1, ..., n}, represents a PSR factor identified through the questionnaire applied and evaluated following the selected instrument. Each PSR factor has a value pj, indicating the attention level required to provide the organization's benefits. The aij cost of attending the j-th PSR factor for the i-th department is defined in the consumption matrix. Finally, ci, i = {1, ..., m} represents the budget limit of each i-th department. According to the optimization model presented in (1) to (3), one backpack represents one company's department. Constraints defined by inequalities (2) indicate that the departments cannot spend more than their budget limits. The objective function (1) is related to maximize the level of attention pj given with the subset of selected PSR represented by the binary variable xj, following the constraints defined by Equation (3). Figure 4 shows the adaptation of the PSR optimization problem as an MKP. Each box represents a PSR using one color, indicating its risk level (see Figure 3). The attention of each PSR generates a cost according to the department where it is applied.
Inside each box, the shapes represent the departments where the risk will be addressed (4 different shapes, one for each department), consuming part of their budget.
The selection of the boxes should generate a maximum benefit of attending to the risks, considering that the sum of the costs of the selected risks should not exceed the budget limit of each department. By using the costs associated with each department, risk levels, and the budget limit assigned to each department in the Case-1 problem, the problem solution is to attend two risks ("Work time management" and "Job Content") in all departments, with a level of attention of 1179 (total profit of v2 and v3). The budget limits are $580, $360, $500, and $380 for the Training, Communication, Industrial Safety, and Human Resources departments. It is interesting to note that if the Training department's budget is reduced to $560, the solution would be to attend to the risk factor "Work time management" only. This change implies reducing the budget use from 96.7% to 53.6%, from 42.4% to 17.6%, from 53.8% to 21.5%, and from 69.7% to 33.0%, in Training, Communication, Industrial Safety, and Human resources departments, respectively.

The SA Algorithm Structure
The SA-based algorithm used for this work is described in the following paragraphs. First, one candidate solution w is defined by the set of selected PSR factors, and the associated costs in each department, as follows: where Profit is the total profit of attending the selected PSRs. xj ∈ {0,1}, for j = {1, …, n}, indicates if a PSR factor is attended or not, and bi, i = {1, …, m}, is the budget consumed by each department attending the selected PSRs.
The cost function f(w) is the total profit of the set of selected PSR factors. The current solution's neighborhood Ω is defined as the set of feasible solutions generated by including or eliminating a selected risk element.
The initial solution, w0, is generated as follows: first, a sequence of n zeros is created. Next, using a random number r ∈ [1, n], the zero value in the r-th position is replaced by a value 1. This replacing process is iteratively applied until the budget limit of any department is exceeded. The resulting binary sequence is used as the initial solution to the SA algorithm.
For constructing a neighboring solution ( Figure 5), a neighborhood function η(w) to swap and shift values on the current solution is defined. Figure 5a shows the swap process, consisting of the generation of two random numbers to obtain a new sequence. Figure 5b shows the shifting process consisting of generating a random number indicating from which element the shifting is performed. A normal random variable with mean 0.5 is used to decide whether swap or shift is applied: if it is greater than 0.5, the swap process is chosen; otherwise, a shift is performed. The SA algorithm adaptation used in this work is shown in Algorithm 1. First, the fSD, Vm, and T0 parameter values are defined (line 1). fSD is a multiple of the standard deviation computed using a set of solutions randomly generated [51]. Vm is a factor to determine the Markov chain size. T0 is the SA initial temperature. T0 is set using the fSD value. Next, the initial solution w0 is created (line 2). Then, in line 3, both current solution w and the optimal solution wopt are set using w0. Additionally, current temperature T is set with T0 and the size of the Markov chain M is computed as follows: where n is the number of PSR factors, and m the company's number of departments. After that, the M value is used to determine the number of iterations of the SA searching procedure. In this procedure, named the Metropolis cycle, new solutions are created and accepted using the Boltzmann criterion (lines [6][7][8][9][10][11][12][13][14][15][16][17][18][19][20]. The current temperature is reduced at the end of each step of the SA external cycle (lines 4-22) using the α coefficient. At the end of the external cycle, the SA algorithm returns its current solution as the best problem solution.
The accepting criterion of new solutions in the Metropolis cycle is as follows: First, a new solution is created (line 7) using the swap and shift operations previously described. Next, the difference ∆ between new and current solutions' fitness values is computed (line 8). Then, if ∆ is greater than 0, the new solution is selected as the current one (line 10). In the other case, if ∆ is not greater than 0, the not improved solution can be selected as the current solution using the Boltzmann criterion (lines 12-15). Generate new solution w' = Ω(w) 8 Calculate if ∆ > 0 then 10 w = w' 11 else 12 Generate random ρ(0,1) Regarding parameter tuning, a sensitivity analysis was performed as follows: for the α parameter, it started from 0.95 to 0.99, where the best was determined to be 0.98. For the β parameter, between 0.01 to 0.001, it was set as 0.001. Two scenarios were tested for the initial SA parameters (T0 and Vm). Table 4 shows the parameter values for these scenarios. The first, called SA-high, T0 = SD and Vm was set between 3 and n, and the second one, named SA-fast, used 0.5 for both parameters. In this way, the solutions generated for all instances are contrasted. A procedure was developed to compute the standard deviation of one hundred random solutions, where the value is recorded in a flat file that is read as input data for the SA algorithm.

Results
In the next paragraphs, the results of the experimental study are detailed. We analyze two test scenarios: first, the results of the HP1 MKP benchmark problem, and then thirty MKP benchmarking instances. This study is carried out on a computer DESKTOP-AQV7GQH, Intel Core i5-8250U CPU, 1.80 GHz, 8.00 GB RAM produced by DELL Technologies, Austin, Texas, USA, using a Windows 10 operating system with a Microsoft Visual C ++ 2012.

Test Scenario for the HP1 MKP Benchmark Problem
According to the MKP model mapped for the PSR optimization problem presented in Section 3, vector p contains the attention levels required by PSR factors, vector x is a binary sequence, where value zero indicates that the PSR is not included, and value one indicates that it is selected. The budget consumption matrix A contains the amounts to attend to each j-th PSR factor of each i-th department, and vector c has the department's maximum budgets. The HP1 instance of the SAC-94 benchmark [54,55] is applied, containing 28 PSR factors with four departments to attend them subject to ceiling budgets. The PSR factors are described in Table 5. A total of 25 PSR factors are used, in accordance with the NOM-035 norm, plus three additional factors from the organizational environment. In Table 6, the values of the risk level for each PSR factor are shown in the third row. Each department's budget limits are shown in the second column of the table, and the first column shows the department names. The budget consumption matrix for each department of the company is shown in the center of Table 6. The amounts corresponding to the V1 factor are presented in column 4 (40,16,38,38), while those for factor V2 are presented in in column 5 (91, 92, 39, 52), and so on. Three solutions to this problem obtained using the SA algorithm are shown in Table  7. Figure 6 shows the percentages of solutions obtained in 30 independent runs. In Table  7, it can be seen that the Human Resources department applies the total of the assigned budget. Table 7. Results for the HP1 benchmark problem using the SA algorithm. The optimal solution is shown in the last row of Table 7, as follows: w = {3418, (1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 1, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1), (216, 199, 201, 180)} In this solution, the attention levels giving the maximum profit of 3418 units selects the following PSR factors: dangerous and unsafe conditions (V1), poor and unsanitary conditions (V2), quantitative loads (V4), accelerated work rates (V5), high responsibility work (V8), lack of control and autonomy over the work (V10), limited or no possibility of development (V11), insufficient participation in management (V12), influence of work outside the workplace (V15), poor clarity of functions (V17), workplace violence (V21), little or no recognition and compensation (V23), limited sense of belonging (V24), job instability (V25), job satisfaction (V26), motivation (V27), and attitude (V28). As shown in Table 6, the departments' budget expenditure is as follows: Training $216, Communication $199, Industrial Safety $201, and Human Resources $180.
As shown in Figure 6, the optimal solution reported for this instance in the existing literature was obtained in 66.6% of the 30 runs executed to the SA algorithm. This proposal dramatically simplifies the preparation of the PSR factors attention plan in a company. When the optimization model is solved, the subset of PSR factors to be addressed is obtained, with a near-optimal benefit level. These factors were adjusted to the budgets of all the departments that correspond to each one. The departments should only generate a budget consumption matrix for each risk factor that they are involved in addressing.
The advantage of using a heuristic such as the SA algorithm to solve the PSR optimization problem is that a set of near-optimal solutions is obtained. Each one is adjusted to the budgets assigned for all departments. These results give the company the possibility of choosing the plan that best suits its needs. In the three solutions shown in Table 7, the solution with 3405 of profit consumes 99% of the budgets of departments 1 and 3, and 98% of that of department 2. This percentage is higher than the optimal solution (3418), but with a budget use of 98% for departments 1 and 2 and 97% for department 3. This fact could be a reasonable justification for the decision maker better implementing the solution with the profit of 3405, indicating a different selection of risk factors only from the following PSR factors: contradictory or inconsistent uploads (V9), insufficient participation in management (V12), poor clarity of functions (V17), and poor relationship with collaborators (V20).

Test Scenario for Thirty MKP Benchmark Problems
Thirty benchmark problems described in Drake et al. [33] and Chu [56] are used to verify the SA algorithm's performance and compare its two sets of parameter values: SAhigh, SA-fast. The SAC-94 library of MKP instances, taken from various articles in the existing literature, represents data from different real-life problems. Most of these instances are small, and their optimal value is known. On the other hand, the GK benchmark [57] presents instances with up to 2500 elements and 100 backpacks. Furthermore, ORLib is a set of 270 instances grouped by the number of elements {100, 250, 500}, number of backpacks {5, 10, 30} and tightness ratio {0.25, 0.50, 0.75}. Most of these studies report only the upper bound obtained. However, in [34,35], global optima were obtained for several OR-Lib instances. These values were obtained using the MIP CPLEX library, but consumed significant amounts of processing time to solve large problem instances. In [34], using a greedy search heuristic considerably reduced the processing time required to find values close to the global optimum, reaching an RE of 0.04%. This heuristic made it possible to discard non-promising search-space areas, using the profit and weights of the MKP elements. The combined use of the greedy approach and Linear Programming (LP) reduced the processing time from 5 h to less than 10 min. The use of metaheuristics for this research makes it possible to solve the problem quickly and easily, with an acceptable RE, allowing its use in real-time decision-making systems in a company. It is also important to clarify that metaheuristics are needed to treat large instances where the processing time grows exponentially, since MKP is classified as an NP-Complete problem [44]. Table 8 shows the characteristics of the 30 benchmarks used in this test scenario. The selected problem size in this work ranges from 28 elements with 4 backpacks up to 250 elements with 30 backpacks. The SA is run 100 times for each of these benchmark instances.  Table 9 shows the main results for both SA-high and SA-fast versions. The first column in this table shows the benchmark instance with the type of PSR optimization problem, where the elements and departments are detailed in Table 8. The second column shows the optimal value for the benchmark instance or the Upper Bound UB* (highlighted using an asterisk) reported in the existing literature. For example, the hp1 instance contains 28 PSR factors attended by four company departments. This table shows the maximum and minimum values reached and the average and mode of the results obtained by running the algorithm 100 times for each instance. According to this table, the mode value is equal to the optimal value in most instances with smaller sizes. It can also be observed that the SA-high version obtained the best results. Alternatively, the SA-fast version equaled the SA-high version results only in the smaller instances and in the OR30X100_0.5_1 instance. Regarding the worst value found (the minimum) by both versions, the SA-fast version found the worst values, except for the following large instances: OR5X250_0.75_4, OR10X250_0.75_4, OR30X250_0.25_3, OR30X250_0.5_3, and OR30X250_0.75_3. With respect to the average values (Avg), the SA-fast version found the worst values for all instances. According to the results shown in Table 9, the SA-high version performance was better than the SA-fast version when maximizing the objective function of the MKP. Table 9. Results for the thirty MKP benchmark problems using the SA algorithm.

Instance
Optimum/ UB* The SA runtime was compared with those obtained in five algorithms used to solve the instances shown in Table 10. The compared were as follows:

SA-High
• CF-LAS (Choice Function-Late acceptance strategy) [33] is a hyper-heuristic with crossover operator performed on an Intel Core 2 Duo 3GHz CPU machine with 2 GB RAM. • RSH (Reduce and Solve Heuristic) [34] coded in C ++ and tested on a PC Intel i5 2.6 GHz CPU, 6 GB RAM. • ICA (Imperialist Competitive Algorithm) [35] also coded in C ++ and run on an AMD Threadripper 2990WX, 3.0 GHz CPU, 64 GB RAM, Windows 10 Pro. • ACO DI (Ant colony optimization with Dynamic Impact) [36] coded in C ++. • 3R-SACRO-PSO (Three-ratio Self-Adaptive Check and Repair Operator-inspired Particle Swarm Optimization) [37] implemented in C and tested using an Intel Core i7 3.4 GHz CPU, 16 GB RAM.
The SA-high and SA-fast versions were implemented in C++, and they were tested on an Intel Core i5 CPU, 8 GB RAM.
The relative error, RE, was evaluated for each instance using Equation (7), where Vop is the optimal global value (the optimum). If this value is unknown, the upper bound reported in the existing literature was used as the reference value. Finally, VSA is the value computed using the proposed SA-based algorithm. RE = abs ((Vop − VSA)/Vop) * 100 (7) Table 10 shows the relative error (RE) obtained by solving 30 MKP instances. The main comparison was performed using the RE, since some of the compared methods do not report their running times. This table shows that the metaheuristics proposed in this work reach competitive results concerning the compared algorithms, highlight the SAhigh version since it obtains the optimal value for the first ten instances. In terms of running times, the SA-fast version was faster than the SA-high version in all instances. Alternatively, comparing the REs between algorithms, the following facts are shown: both SA versions were better than the CF-LAS method when solving the instance number 10. Furthermore, in 18 large instances (instances [13][14][15][16][17][18][19][20][21][22][23][24][25][26][27][28][29][30], the SA-high version was better than the CF-LAS method. In particular, the SA-high performance was better than RSH in one instance and obtained the same values in two instances. In the rest of the cases, both SA versions were competitive with all metaheuristics. The highest RE obtained for the SAhigh version was less than 1, and for the SA-fast version, this was at least 1.12. These results indicate that the SA versions proposed in this work have outstanding efficiency, and they compete with the metaheuristics found in the literature for the MKP problem. In Table 10, column 2, the optimum is the global optimum in most cases, and in very few cases, the optimum is an upper bound (marked with *), because the global optimum of the MKP is not yet known.  Figure 7 shows the relative error of the thirty benchmark instances. The highest value for the SA-high scheme is 0.93 and that for the SA-fast scheme is 1.12. For benchmark instance number 10, the SA-high scheme reached the global optimum, while the SA-fast scheme had a relative error of 0.52. Carrying out one contrast between the two SA parameter schemes, Figure 8 shows that, on average, the SA-fast scheme solved the benchmark instances 76% faster than the SA-high scheme.

Statistical Analysis
The algorithms CF-LAS, RSH, and the SA-high version were used to perform the statistical analysis of the algorithm proposed in this work. The comparison using these algorithms was performed on the basis of Table 10, since the others were lacking in data. For the statistical analysis, the results of 18 large instances were used, from the OR5X100_0.25_6 instance to the OR30X250_0.75_3 instance (Table 11). The null hypothesis H0 used in this analysis indicates that the means of the results of the algorithms are equal, and the alternative hypothesis H1 suggests that at least one algorithm had a different behavior.
:  Figure 9 presents the normality graph of the three algorithms. From this figure, it can be concluded that no normality exists in the data since the points are not located on the diagonal of the graphs. The homoscedasticity analysis was performed using the box-and-whisker graph presented in Figure 10. It can be seen that, for each algorithm, the boxes are not the same, indicating a difference in variances and that homoscedasticity does not exist. As the two assumptions of normality and homoscedasticity of the data did not exist for applying a parametric test such as ANOVA, the robust ANOVA test of Welch and Box [58] was used. The Welch and Box test [58] was performed in order to perform statistical analysis with robust ANOVA using the WRS R package [59] to compare the three algorithms using the 10% of the trimmed means. The obtained p-value was 1.912572 × 10 −5 , and the null hypothesis was refuted. The value Fw of the statistician was 18.31316, and the degrees of freedom were v1 = 2 and v2 = 22.47213.

Discussion
This methodology is used to generate the first plan of attention to PSR factors. In the case study, the budget consumption matrix is obtained from the hp1 instance. Still, in a real case, its preparation requires great effort, since multiple meetings between the company departments are needed to agree on the amounts used to attend each risk factor and establish the budget ceiling. When solving the optimization model of this instance, Table  7 shows the three solutions obtained using the budget consumed by each department. Again, it is essential to emphasize that companies must apply most of their assigned budget in order to achieve the most significant benefit, since an unconsumed budget would imply neglecting its purpose or mismanagement. In the instance of hp1, a competitive solution achieved a profit of 3405, consuming 99% of the budget of departments 1 and 3 and 98% of the budget of department 2. The benefits of using metaheuristics such as the SA-based approach presented in this work are evident. They can obtain several near-optimal solutions, and the decision-makers have several scenarios from which to select that which best aligns with the company's objectives, such as the use of the greatest budget to obtain the greatest benefit in addressing psychosocial risks.
The cost matrix generation for each department, even without considering the risk levels presented in the company, provides the opportunity to incorporate future actions and move towards the continuous improvement of the organizational environment. In accordance with various methodologies described in the existing literature, the success of these projects type (continuous improvement) lies in communication. However, since this work involves psychosocial aspects, the application of diverse issues such as social support, the use of mechanisms of disagreement, and well-targeted training will make a difference. This methodology requires costs to attend to each PSR factor offering important support in decision making. This proposal reduces the tensions related to the creation of one intervention plan, since each department generally has its own prevention perspective according to its own goals and objectives. The generation of the cost matrix per department enriches the actions, activities and interventions that need to be followed in order to improve the organizational environment.
Regarding the metaheuristic used to solve the PSR optimization problem, it is clear that the smallest instances, the solution found by the SA algorithm, is the optimal value reported in the benchmark literature. For the other instances, the sampling error shown in Figure 8 validates the two SA schemes' quality and precision. The SA-fast scheme has a maximum RE no greater than 1.12%, representing competitive solutions found in a reasonable time. It is ideal for large PSR optimization problems. This behavior allows the SA algorithm to be incorporated into the technology platform with the enterprise security methodology [60][61][62]. On the other hand, as shown for small instances, the SA-high scheme presents greater precision, where its implementation brings great benefits.
Regarding the statistical results, we observe that the proposed SA algorithm is competitive with those contrasted in the existing literature, demonstrating its independence.
Solving the PSR optimization problem has great relevance, due to the changes in the forms of work that are currently taking place [1,[7][8][9]63]. Algorithmic management [2,[64][65][66], which assigns tasks and evaluates workers using data-driven systems, has proven to be efficient. However, at the same time, the present challenges of the psychosocial type should be considered. This presents a wide field of study for future works.
As future work, we intend to automate the identification of the level of attention that the PSR factor requires, incorporating, in addition to the questionnaire, information from interviews, minutes of discussion group meetings, and definitions of performance indicators aligned to the strategic objectives of the company, among other things.

Conclusions
Given the challenges that organizations face regarding digitization and process automation, technological tools play an essential role in supporting decision making. Changes in the workforce go hand in hand with job wellness and occupational health, hence the importance of using algorithms that optimize solutions to address psychosocial risks. In this work, a methodological scheme was presented that, on the basis of the detection of psychosocial risk factors in a company, the mapping to an MKP optimization model, and its solution using the SA algorithm, is able to obtain the subset of risk factors for a company with a maximum value at the level of care, that is, with a local-optimal solution that maximizes the benefit to the factors with the highest identified risk. This subset of selected factors can be incorporated into an optimal PSR treatment plan in the workplace. As shown by the evaluation of the method, the level of attention to risks is maximized under the limited budgets of a company's departments. This methodology favors the organizational environment and promotes business competitiveness, complying with existing regulations for the identification, analysis and prevention of PSR factors.
By solving the benchmark instances for MKP, it was shown that the SA algorithm was able to solve the PSR optimization problem. The algorithm obtained values close to the optimum (and, in some cases, the global optimum) in the conducted tests.
The suggested metaheuristic was developed using two tuning schemes: SA-high and SA-fast. It can be observed that the second scheme had a better execution time (almost 80% faster than the first), sacrificing precision of results by a value within the range of 2-3%. These results are competitive, since the operation speed is one of the variables of interest under the new requirements given in the work scheme, such as algorithm management.
The results show that an SA-high scheme is a promising approach, since it achieves results below 1% for relative error. This implies that the results may be beneficial for a real PSR problem of a company.
As future work, it is intended to automate the identification of the level of attention that the PSR factor requires, incorporating, in addition to the questionnaire, information from interviews, minutes of discussion group meetings, and the definition of performance indicators aligned to the strategic objectives of the company, among other things.