To study the algorithm, I tried to solve the shift scheduling problem. The hottest topic in May 2020, when I was writing this article, was COVID-19, and when I think about how this will affect my company's livelihood, I would like to come to work (requests from employees. Until last year. I think it's a good idea to satisfy both the company's request to go to work) and the telework (company's request. Until last year, the employee's request to telework). Yeah, it's a so-called shift scheduling problem.
First, let's try the Ising model, which is also used in D-Wave's quantum annealing method. However, quantum annealing cannot be executed on an ordinary computer, so we will use ordinary annealing. It's the same way as Fujitsu's Digital Annealer.
However, it was troublesome to create a process to solve the Ising model by simulated annealing, so [D-Wave's neal](https://docs.ocean.dwavesys.com/projects/neal/en/latest/index. I used html). Also, it is quite difficult to make the Ising model by hand (or rather, my math ability is too low to understand the fragments), so the generation of the Ising model was [PyQUBO of Recruit Communications Co., Ltd.](https: // I used pyqubo.readthedocs.io/en/latest/).
So, the code to solve the problem using neal and PyQUBO looks like this.
import numpy as np
from funcy import *
from neal import SimulatedAnnealingSampler
from pyqubo import Array, Constraint, Placeholder, Sum
M = 5 # Number of employees
D = 10 # days
BETA_RANGE = (5, 100) # Reciprocal of annealing temperature. The larger the solution, the more stable the solution, but the more likely it is to fall into a local solution.
NUM_READS = 10 # Number of annealings. NUM_READS solutions are generated. The more you have, the more likely you are to get a better solution.
NUM_SWEEPS = 100000 # The number of times to perform the annealing steps. The number of iterations to generate a single solution. The larger the value, the more likely it is that a better solution will be obtained.
BETA_SCHEDULE_TYPE ='linear' # How to change the annealing temperature. If it is linear, it will change linearly.
# Anneal the Ising model using neal and return the solution
def solve(hs, js):
response = SimulatedAnnealingSampler().sample_ising(hs, js, beta_range=BETA_RANGE, num_reads=NUM_READS, num_sweeps=NUM_SWEEPS, beta_schedule_type=BETA_SCHEDULE_TYPE, seed=1)
# NUM_READS Returns the best solution out of
return tuple(response.record.sample[np.argmin(response.record.energy)])
# Define the variables that make up QUBO
xs = Array.create('x', shape=(M, D), vartype='BINARY')
# Define variables for tuning
a = Placeholder('A')
b = Placeholder('B')
# Define QUBO. from here……
h = 0
# Add a constraint of 2 or more people per day and as few as possible. Penalties are now incurred for more or less than 2 people
for d in range(D):
h + = a * Constraint ((Sum (0, M, lambda m: xs [m] [d]) --2) ** 2, f'day- {d}') # 2 is negative if at least , If there are many, it will be a positive number, but square it to make it a positive value
# Add a constraint that you will not come to work with the same person on another day
for m1 in range(M):
for m2 in range(m1 + 1, M):
for d1 in range(D):
for d2 in range(d1 + 1, D):
h + = b * xs [m1] [d1] * xs [m2] [d1] * xs [m1] [d2] * xs [m2] [d2] # xs is 1 or 0, so if you want to multiply It will be 1 only if all are 1.
# Compile and make a model
model = h.compile()
# ……So far. Define QUBO
# Variable values for tuning
feed_dict = {'A': 2.0, 'B': 1.0}
# Generate an Ising model and solve it with neal
hs, js, _ = model.to_ising(feed_dict=feed_dict)
answer, broken, energy = model.decode_solution(dict(enumerate(solve(hs, js))), vartype='SPIN', feed_dict=feed_dict)
# Output the result
print (f'broken: \ t {len (broken)}') # If Constraint is violated, broken will be populated.
print (f'energy: \ t {energy}') # QUBO's energy. In this model, it will be 0 if all constraints are met.
# Outputs employees who come to work on a daily basis
for d in range(D):
print(tuple(keep(lambda m: 'ABCDE'[m] if answer['x'][m][d] == 1 else False, range(M))))
So, the result of executing this program is like this.
broken: 0
energy: 0.0
('D', 'E')
('A', 'D')
('B', 'C')
('A', 'C')
('B', 'D')
('B', 'E')
('C', 'E')
('A', 'E')
('A', 'B')
('C', 'D')
Yeah, that's right. Two people come to work every day, and there is no same combination. By the way, the execution time on my computer was 0.734 seconds. Since 10 solutions have been created, 0.073 seconds per solution. neal is really fast!
[Wikipedia Ising model page](https://ja.wikipedia.org/wiki/%E3%82%A4%E3%82%B8%E3%83%B3%E3%82%B0%E6%A8% If you read A1% E5% 9E% 8B), you can understand that it is a lattice model that is composed of lattice points that take two coordination states and considers the interaction of only the adjacent lattice points. not. Hamiltonian is a physical term, but is it delicious? Do I have to study physics from now on?
Please be assured. It is the people who make the hardware that utilizes the Ising model who have to study physics, not us programmers. Moreover, the Ising model is a fairly simple story in terms of software.
Let me explain with a concrete example. First, ignore the word grid (which seems to be vertical and horizontal in the image) (important for physics hardware stores, but not for us programmers). So, think of the two coordination states as special variables that can only be assigned either -1 or 1. It's hard to get an image of 1 or -1, so let's consider the case where three people vote for or against a proposal. If you agree, it will be 1, and if you disagree, it will be -1. There is a delicate relationship between these three people, and it feels good if Mr. a and Mr. b take the same action, and it feels good if Mr. c and Mr. b take different actions due to some past cause. I will. Now, how can we make everyone feel as comfortable as possible under these conditions?
Computers can only handle numbers, so somehow they have to be expressed numerically. Since a, b, and c can only take either 1 or -1, expressions that make the numbers * smaller * when everyone feels comfortable can be expressed with this code.
-1 * a * b + 1 * b * c
Let's take a closer look. The result of multiplying the combination of 1 and -1 is as follows.
| Value 1 | Value 2 | result |
|---|---|---|
| 1 | 1 | 1 |
| 1 | -1 | -1 |
| -1 | 1 | -1 |
| -1 | -1 | 1 |
When the values are the same and different, it is comfortably divided into 1 and -1, so if you multiply this result by -1, it will be a small value of -1 if it is the same, and a large value of 1 if it is different. You see, you could express it well with the above formula, right?
However, it is inconvenient because it is not versatile as it is, so let's make it versatile. Add the values in the table below to the multiplication results of all variables (however, a * a is the same every time, so it is meaningless, so leave it blank. I will hang it.
| a | b | c | |
| a | -1 | 0 | |
| b | 1 | ||
| c | |||
You can set it to 0 where it doesn't matter. So, when I write the code to calculate using this table, it looks like this. Imagine that a, b, and c are in an array called xs, and the values in the table above are in a variable called js.
def energy():
result = 0
for i in range(1, len(xs)):
for j in range(i + 1, len(xs)):
result += js[i][j] * xs[i] * xs[j]
return result
Now, even if the relationship changes due to reconciliation between Mr. b and Mr. c, it can be expressed with the same code (actually, to make it more general (ʻaand). We also use another variable (hs in Simulated Annealing code using the Ising model) (to take advantage of whether b or c` is 1 or -1). And this versatile and wonderful is the Ising model.
At this point, it's easy to see if it's better or worse by flipping the appropriate element of xs and comparing the results in the code above. It can be solved by the so-called mountain climbing method. However, when I looked it up, the explanation of the hill climbing method said that it was easy to fall into a local solution, so I was a little worried, so even if the result of the above code is a little worse at the beginning, I will allow movement, the last one is that degree Let's work around this problem by reducing and moving it in the direction of improvement. This is called Simulated Annealing because it resembles Simulated Annealing in the real world, and the degree to which it allows even if the result is poor is called temperature, following Simulated Annealing in the real world. Also, since it is energy that corresponds to the result of the above code in the real world, the result of the above code is called energy. Therefore, it is an expression like finding a combination that minimizes energy while lowering the temperature. D-Wave's neal does this very efficiently.
So if you anneal it, the result should be a = 1, b = 1, c = -1 or a = -1, b = -1, c = 1. In both cases, the energy is the minimum -2, so everyone feels good.
Then all you have to do is create the values in the table above according to the constraints mentioned at the beginning of this article ... but this is quite difficult. I don't know what to do when multiple variables are involved, as in the constraints of this article (it seems that you add a variable behind the scenes and define it in relation to the added variable). So let's use Recruit Communications' PyQUBO. With PyQUBO, the expression in the formula is converted into a versatile Ising model. Also, since it is difficult to think of 1 and -1, it can be defined in the form of QUBO using 1 and 0 (I could not understand, but it seems that the Ising model and QUBO can be converted to each other). Ising.
Let's express it concretely. If you think of xs as a two-dimensional array of M × D and set it to 1 if you go to work and 0 if you do not go to work, the number of employees who come to work on the first day can be calculated as follows.
result = 0
for m in range(M):
result += xs[m][0]
return result
So, the closer this result is to 2, the better, and it doesn't matter if it's more or less. So I want to subtract 2 and do ʻabs to find the absolute value ... but unfortunately I can't use a function like ʻabs. So what to do is to square it. I'm not good at math, so I'm not sure, but it seems that multiplying minus by minus will give you a plus. Also, PyQUBO provides a function to calculate the sum called Sum, so the code is as follows.
h += (Sum(0, M, lambda m: xs[m][0]) - 2) ** 2
Now, when the number of people who come to work is 2, the energy is the smallest, and the energy is greater than or less than 2.
For the rest of the different pairs of employees, it's okay if you don't have the same combination on different days. For the sake of simplicity, if you think in a little more detail, for example, "Employee 0 and Employee 1 will be penalized if the 2nd and 3rd days are the same", you can express it with the following code.
h += xs[0][2] * xs[1][2] * xs[0][3] * xs[1][3]
Since it is QUBO, the value of the element of xs is either 1 or 0, so it will be 0 unless both employee 0 and employee 1 come to work on both day 0 and day 1, that is, all 1's. .. Since all 1s are the same combination, it is very convenient as a penalty if it becomes 1 in that case. All you have to do is write a quadruple loop to cover the combinations of other employees and the combinations of other days.
Also, most things have priorities, and it's difficult to add different things together (I don't think the unit systems of the two well-made formulas this time are the same: "itch" and "it". Isn't it difficult to calculate "unpleasantness" by adding "noisiness"?). So, let's multiply the appropriate coefficients, ʻaandb, so that these values can be specified at the time of tuning (itchiness x ʻa + noisy x b = unpleasantness and provisional. Then, adjust the values of ʻaandbas appropriate).Placeholder is convenient in such a case. This time, I tried various things and set ʻa = 1.0 and b = 0.5. If you define QUBO with PyQUBO like this, it will be converted to Ising model in one shot with model.to_ising (), and if you solve it with neal, the answer will be returned. Also, you can use D-Wave instead of neal. You can also solve it with a digital annealer and compare the execution time and accuracy with neal. Well, it has become convenient in the world.
Let's get on with it and do it with a genetic algorithm that makes the name feel romantic.
As usual, it was tedious to create a process to do a genetic algorithm, so I used an open source library called DEAP. The code looks like this.
from deap import algorithms, base, creator, tools
from functools import reduce
from funcy import *
from random import randint, random, seed
M = 5 # Number of employees
D = 10 # days
# Evaluation function
def evaluate(individual):
# Add a constraint of 2 or more people per day and as few as possible. Penalties are now incurred for more or less than 2 people
def member_size():
result = 0
for d in range(D):
result + = abs (reduce (lambda acc, m: acc + individual [m * D + d], range (M), 0) --2) # Since the value itself is used, abs can also be used.
return result
# Add a constraint that you will not come to work with the same person on another day
def different_member():
result = 0
for m1 in range(M):
for m2 in range(m1 + 1, M):
for d1 in range(D):
for d2 in range(d1 + 1, D):
result += individual[m1 * D + d1] * individual[m2 * D + d1] * individual[m1 * D + d2] * individual[m2 * D + d2]
return result
# Returns multiple evaluation viewpoints as a tuple with the evaluation results from each viewpoint as an element.
return (member_size(), different_member())
# DEAP defines how to genetically algorithm
creator.create ('Fitness', base.Fitness, weights = (-1.0, -0.5)) The smaller the result of # evaluate (), the better, so add a minus to the weight.
creator.create('Individual', list, fitness=creator.Fitness)
toolbox = base.Toolbox()
toolbox.register('attribute', randint, 0, 1)
toolbox.register('individual', tools.initRepeat, creator.Individual, toolbox.attribute, n=M * D)
toolbox.register('population', tools.initRepeat, list, toolbox.individual)
toolbox.register('mate', tools.cxTwoPoint)
toolbox.register('mutate', tools.mutFlipBit, indpb=0.05)
toolbox.register('select', tools.selTournament, tournsize=3)
toolbox.register('evaluate', evaluate)
# Random seeds are fixed for reproducibility
seed(1)
# Solve the problem with a genetic algorithm
population, _ = algorithms.eaSimple(toolbox.population(n=100), toolbox, 0.5, 0.2, 300, verbose=False)
# Get the best solution
individual = tools.selBest(population, 1)[0]
# Output the result
print(f'fitness:\t{individual.fitness.values}')
# Outputs employees who come to work on a daily basis
for d in range(D):
print(tuple(keep(lambda m: 'ABCDE'[m] if individual[m * D + d] else False, range(M))))
And the result is like this.
fitness: (0.0, 0.0)
('B', 'E')
('A', 'D')
('B', 'C')
('C', 'E')
('D', 'E')
('C', 'D')
('A', 'C')
('A', 'E')
('B', 'D')
('A', 'B')
Yes. That's the correct solution. The execution time on my computer was 4.595 seconds. If you feel that it is too slow to use, please read the following "Simple explanation" to the end.
A child born of a cool parent is very cool. And I was born from a parent who wasn't ...
As with the previous Simulated Annealing method, it is important to create a new solution in the method of searching for a better solution by creating a new solution. For example, if you randomly create a new solution, it will probably get worse and you will not find a solution that you can wait for. Therefore, the hill climbing method uses a solution in the vicinity of the current solution, which is a slight modification of the current solution. I think it's probably good because it resembles a good solution. So, in the genetic algorithm, the solution is expressed like a gene, a new solution is created by mating and a child is born or mutated, and a better solution is created by natural selection. I will look for it. For mating and natural selection, we have multiple solutions instead of one. The children of cool parents are pretty cool, so I feel like I can survive without being eliminated in the marriage market.
So, the most important thing in the genetic algorithm is "how to express the solution with chromosomes". Like this time, if you come to work and do not come to work with 1, you can express it with a list of 0. Any set of numbers can be used (NumPy's ndarray, Set, Dictionary, tree structure, etc. can be used), so for example, if it is a traveling salesman problem, a list of the numbers of the cities to be visited is okay. When driving a car, the strength with which the accelerator or brake is depressed may be expressed in floating point. DEAP has a wealth of functions that enable the expression of various genes and chromosomes. For example, if the number of the city to be visited by the traveling salesman is a gene, Creating Types of DEAP documentation Permutation is useful (now you don't have to use tedious techniques like ordinal representation anymore!).
Also, there are actually various methods such as mating (called crossing in the genetic algorithm), mutation, and natural selection, but they implement many of them. And it also provides a nice way to combine these in the ʻalgorithms` package.
However, DEAP has a little quirk in how to use it ... We will reuse the functions provided by DEAP to create the tools necessary to solve the problem, but we have to do it with the functions of DEAP instead of ordinary function synthesis. For example, to define an individual (corresponding to one of the solutions) ʻIndividual class, you have to do it with the DEAP API like creator.create ('Invididual', ...) . Hmm. So, creating intersecting methods is like toolbox.register ('mate', tools.cxTwoPoint). It's easy to generate a method that intersects two points with this alone, but if possible I wanted to write it like a normal partial` ...
Well, this is a luxury problem, so let's make a program quickly. Since you have to write the function to evaluate the individual by yourself, you can refer to the code written at the time of annealing using the Ising model, but in the case of the genetic algorithm, ʻabs can be used and it is convenient. While thinking, I created the ʻevaluate () function. After that, create a Fitness class that evaluates the goodness of the solution, a ʻIndividual class that expresses an individual, create a ʻattribute () method that creates the attributes of the individual's chromosome, and use it to create an individual. Create a ʻindividual ()method that creates apopulation ()method that uses it to generate a population. After that, the intersectionmate ()method required for the genetic algorithm, the mutationmutation ()method, the natural selectionselect () method, and the first created ʻevalute () function Generate a ʻevalute ()` method to call.
So, this time, I tried to execute the simplest form of genetic algorithm with ʻalgorithms.eaSimple ()`. Even if you leave it to the library completely, if you do a genetic algorithm for 300 generations with 100 individuals, you will get the correct answer.
By the way, the good points of the genetic algorithm, which was slower than Simulated Annealing using the Ising model, are that there are many applicable problems because the expression of the solution is more flexible than the Ising model, and there is a lot of room for tuning because there are many methods. is. I didn't tune it in this article, but if you make the chromosome representation more efficient, change the way it crosses, or change the probability of mutations to be nice, it will be more than simulated annealing using the Ising model. It may be possible to derive a highly accurate solution at high speed. Tuning is fun because you can see the effect immediately.
Well, in order to do that, I have to study various methods of genetic algorithms, but if I study while actually trying with DEAP, I think I can master it immediately.
When I wondered if there was anything else, I noticed the integer programming method. It's a version that makes Linear Programming even more difficult by limiting it to integers.
Linear programming is a linear expression (an expression in which only one variable is multiplied and connected by addition and subtraction. 3 * x + y is a linear expression, and 3 * a * a + b seems to be a quadratic expression). It is a method to express an objective function or a constraint function and solve it mathematically quickly. The objective function is an expression that expresses the goodness of a solution as if it was done by simulated annealing using the Ising model or a genetic algorithm, and selects a solution that is as large or small as possible. So, the constraint function is a conditional expression that expresses the constraint of the solution.
e? Do you know what you are saying?
I don't know at all so don't ask ... However, if you use PuLP, you can do integer programming (and of course linear programming) quickly. The code looks like this.
from functools import reduce
from funcy import *
from pulp import *
M = 5 # Number of employees
D = 10 # days
# Define the variables used in the problem
xs = LpVariable.dicts('x', (range(M), range(D)), 0, 1, 'Binary')
# Define the problem. from here……
problem = LpProblem('shift-scheduling-problem', LpMinimize)
# Add a constraint of 2 or more people per day
for d in range(D):
problem += reduce(lambda acc, m: acc + xs[m][d], range(M), 0) >= 2
# Add a constraint that you will not come to work with the same person on another day
for m1 in range(M):
for m2 in range(m1 + 1, M):
for d1 in range(D):
for d2 in range(d1 + 1, D):
problem + = xs [m1] [d1] + xs [m2] [d1] + xs [m1] [d2] + xs [m2] [d2] <= 3 # The inequality sign (<) could not be used, so < = 3
problem.writeLP('shift-scheduling-problem')
# ……So far. Define the problem
# Solve the problem with integer programming
status = problem.solve()
# Output the result
print(LpStatus[status])
# Outputs employees who come to work on a daily basis
for d in range(D):
print(tuple(keep(lambda m: 'ABCDE'[m] if xs[m][d].value() else False, range(M))))
And the result is like this.
Optimal
('D', 'E')
('C', 'D')
('A', 'E')
('C', 'E')
('B', 'D')
('B', 'C')
('A', 'C')
('A', 'D')
('B', 'E')
('A', 'B')
It's an optimal solution. The execution time is 0.132 seconds, which is very fast!
Variables in PuLP are created with LpVariable. This time I needed a lot of variables, so I generated a lot with LpVariable.dicts () at once. Also, this time, there is no superiority or inferiority in the solution that satisfies the constraint (if you try to satisfy the constraint of different pairs of employees as much as possible, the number of people who come to work will decrease), so only the constraint is defined.
The way to write constraints in PuLP is a conditional expression. 1 The result of calculating the number of people at Hinode with reduce () is like > = 2. We will add this conditional expression to the problem created as LpProblem with+ =. With the constraint that you do not come to work with the same person on another day, if you do xs [] * xs [] * xs [] * xs [] as before, it will be a quaternary expression, so add (in this case) As a result of (linear expression), I made a constraint in the form of <= 3.
All you have to do is solve () this. If the constraint can be satisfied, the optimal solution with the smallest objective function in the constraint can be solved by mathematical magic and returned (By the way, you can check whether the constraint is satisfied by LpStatus [status]. Masu). I don't know what kind of math magic I'm using, but it's okay if I just "use" PuLP.
Yeah, I don't know what's inside, but PuLP is amazing to find the optimal solution in such a short time ... Of course, it's amazing, but unfortunately, integer programming is not perfect. If it's as easy as this one, you'll get an answer right away, but if it's a difficult one, it can take a very long time to find a solution. Simulated annealing using a genetic algorithm or Ising model is a method that outputs a solution that seems to be reasonably good, although it may not be optimal, so even difficult problems can be solved in some way. If it is a difficult problem, it may not be possible to express it with a linear expression. Of course, if the solution is not absolutely optimal, or if the problem is not very complicated, the integer programming method (linear programming method) is better.
But that? If you calm down, you might be able to solve it in a shorter time with a simpler program. See, thinking in combination, it looks like this …….
from itertools import combinations, cycle
# Generate a unique combination of two employees
members = cycle (combinations ('ABCDE', 2)) # You don't have to cycle for 10 days, but just in case
# Output for 10 days
for _ in range(10):
print(next(members))
If you try it ...
('A', 'B')
('A', 'C')
('A', 'D')
('A', 'E')
('B', 'C')
('B', 'D')
('B', 'E')
('C', 'D')
('C', 'E')
('D', 'E')
Yeah, it's correct to see. The execution time was 0.028 seconds ....
In fact, the process of choosing different pairs is very easy if you can express it programmatically. Like this.
def getPairs(xs):
for i in range(len(xs)):
for j in range(i + 1, len(xs)):
yield xs[i], xs[j]
I just started the range of the j loop with ʻi + 1` so that the same ones aren't chosen, and the combinations in reverse order aren't chosen. It's the same technique used in Simulated Annealing code using the Ising model, where the same person doesn't come to work on different days. So, when I counted the results returned by this function, it was 10, so the problem this time was that if done well, I could get an answer that meets all the constraints.
So, the number of combinations to choose 2 out of 5 is 10, I feel like I learned somewhere in the past ... When I looked it up, [Wikipedia Combinatorics](https://ja.wikipedia.org/wiki/%E7%B5%84%E5%90%88%E3%81%9B%E6%95%B0%E5 That is exactly the combination formula that does not allow repetition of% AD% A6). 5! / (2! * (5-2)!) = 10.
And, because it is a famous process like this, the combinations are libraryed in most programming languages. For Python used in this article, that is ʻitertools.combinarions. It was okay to just apply the result of combinations to cycle`, which repeats the set to create an infinite set, and display the number of days from the beginning.
So this is the punch line this time (I'm sorry for those who noticed the punch line due to constraints. Also, the term shift scheduling problem is fishing. I'm really sorry for those who are seriously doing shift scheduling problems). But what I want to say in this article is not about mastering combinations. This paper argues that simulated annealing using the Ising model, genetic algorithms, integer programming, and combinations can be easily implemented using modern libraries. Each method has its advantages and disadvantages, so the problem depends on which method is appropriate. Then, when asked what to do, I think I should try various things for the time being. Because it's so easy to implement.
Recommended Posts