遺伝的アルゴリズムによるスケジューリング問題の解法

こんにちは，機械学習エンジニアの大下です．
就活がなかなか上手くいかなくて，つらいですね．
お仕事ください．

今回は，スケジューリング問題を扱いたいと思います．スケジューリング問題には例えば，
看護師のシフトを組む，ナーススケジューリング問題．それらを一般化した○○（バス，鉄道，航空など）勤務スケジューリング問題があります．

##スケジューリング問題とは
スケジューリング問題とは最適化問題の一種です．最適化問題には，例えば（非）線形計画問題，整数計画問題，半正定値計画問題，離散最適化問題（組み合わせ最適化問題）などの問題があります．
今回は最適化問題の中のスケジューリング問題に限って説明します．

##進化計算アルゴリズムとは
進化計算アルゴリズムとは遺伝的アルゴリズム、遺伝的プログラミング、進化戦略、進化的プログラミングを総称した呼び方です。

今回はそのなかでも遺伝的アルゴリズムを使用するので、その説明をしたいと思います。

##遺伝的アルゴリズムとは
遺伝的アルゴリズム(Genetic Algorithm:GA)とは1960年代にHollandによって提案された生物の進化過程を模擬するアルゴリズムである．
簡単に生物の進化過程を説明する．
生物は親の情報を部分的に受け継がれる．その過程には細胞が分裂する際に細胞の中にある遺伝子(gene)が複製されることによって実現される．細胞分裂時には染色体(chromosome)として多数の遺伝子が棒状にまとまる．各遺伝子の染色体ないでの位置を遺伝子座(locus)と呼ぶ．また遺伝子には複数の種類が存在し，それを対立遺伝子（allele）と呼ぶ．これらの対立遺伝子の並びを遺伝子型(genotype)と呼び，遺伝子型によって決定される形質を表現型(phenotype)と呼ぶ．

そして遺伝子の位置を交換する交差，一部の遺伝子を変化させる突然変異，などによって進化過程をモデル化しているのが遺伝的アルゴリズムである．

 

##ナーススケジューリング問題とは
ナーススケジューリング問題(Nurse Scheduling Probrem:NSP)とは，看護師の勤務日を制約条件を満たすように勤務日の配置を考える問題である．

##制約条件とは
制約条件とは不等式の形として与えられる．
例えば看護師の勤務時間をTとすると

```python
if T >= 8:
  penalty = penalty + 1
```

のように実装する．（コード参照）

##適応度とは
その個体がどれほど環境に適応しているかを表す値である．
スケジューリング問題ではペナルティ関数の値と考える．

###適応度地形とは
適応度地形（英語: fitness landscape）は1923年にシューアル・ライトらによって提唱された、遺伝子型と生殖成功率の関係（適応度）を視覚化するために用いられる数学的モデルである。(wikipedia)
<blockquote class="twitter-tweet" data-lang="ja"><p lang="ja" dir="ltr">Using fitness landscapes to visualize evolution in action <a href="https://t.co/pDG9gKsfro">https://t.co/pDG9gKsfro</a> <a href="https://twitter.com/YouTube?ref_src=twsrc%5Etfw">@YouTube</a>さんから</p>&mdash; NoriakiOshita (@whisponchan) <a href="https://twitter.com/whisponchan/status/1137785936631275530?ref_src=twsrc%5Etfw">2019年6月9日</a></blockquote>
<script async src="https://platform.twitter.com/widgets.js" charset="utf-8"></script>

###実際の問題に適用する

オリジナルの勤務表データを作成して、それについてGAで最適化を行う

(次回以降にします．)

この問題では勤務時間帯は朝，昼，夜，深夜の4パターンのみとします．
そして7人分の一週間のシフトを自動作成するという限定的な問題に限ります．
雰囲気だけを掴んで貰うためこのような条件にしています．


![qiitaのGAの記事.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/33329/388373f8-13d0-b1ae-e8fe-60d821266ac8.jpeg)

###解の評価
ペナルティ関数を最小することによって解は局所最適解に求まります．


##ソースコード
ソースコードはこちらのサイトを参考にしました．
https://qiita.com/Azunyan1111/items/975c67129d99de33dc21

```main.py
    class Genom:

      genom_list = None
      evaluation = None

      def __init__(self, genom_list, evaluation):
          self.genom_list = genom_list
          self.evaluation = evaluation


      def getGenom(self):
          return self.genom_list


      def getEvaluation(self):
          return self.evaluation


      def setGenom(self, genom_list):
          self.genom_list = genom_list


      def setEvaluation(self, evaluation):
          self.evaluation = evaluation

# -*- coding: utf-8 -*-
import random
from decimal import Decimal
import numpy as np
from itertools import zip_longest

# シフト表の日数
DAY = 7
# 遺伝子情報の長さ(4*people*day)
GENOM_LENGTH = 4 * 30 * DAY
# 遺伝子集団の大きさ
MAX_GENOM_LIST = 100
# 遺伝子選択数
SELECT_GENOM = 20
# 個体突然変異確率
INDIVIDUAL_MUTATION = 0.1
# 遺伝子突然変異確率
GENOM_MUTATION = 0.1
# 繰り返す世代数
MAX_GENERATION = 150

# 朝、昼、夜、深夜帯の規定された時間範囲、単位は時間
element = [3.0, 2.0, 3.0, 8.0]
diagram = []
[diagram.append(element) for m in range(int(GENOM_LENGTH / 4 / DAY))]
diagram = np.reshape(diagram, (int(GENOM_LENGTH / 4 / DAY), 4))

def create_genom():
    """
    引数で指定された桁のランダムな遺伝子情報を生成、格納したgenomClassで返します。
    :param length: 遺伝子情報の長さ
    :return: 生成した個体集団genomClass
    """
    genome_element_list = []
    genome_days_list = []
    genome_all_list = []

    for i in range(int(GENOM_LENGTH / 4 / DAY)):
        for j in range(DAY):
            for k in range(4):
                genome_element_list.append(float(random.randint(0, 1)))
            genome_days_list.append(genome_element_list)
            genome_element_list = []
        genome_all_list.append(genome_days_list)
        genome_days_list = []
    return Genom(genome_all_list, 0.0)

penalty = 0.0

def evaluation(ga):
    """評価関数です。今回は全ての遺伝子が1となれば最適解となるので、
    合計して遺伝子と同じ長さの数となった場合を1として0.00~1.00で評価します
    :param ga: 評価を行うgenomClass
    :return: 評価処理をしたgenomClassを返す
    """
    global penalty
    penalty = 0.0
    # 多次元配列から行列にする
    ga_array = ga.getGenom()

    penalty = evaluation_function(ga_array)
    return penalty

def evaluation_function(ga_array):
    global penalty
    sum_time = 0.0
    for i in range(len(ga_array)):
        for j in range(len(ga_array[i])):
            sum_time = 0.0
            for k in range(4):
                #print(ga_array[i][j])
                if type(ga_array[i][j]) == list:
                    if ga_array[i][j][k] == 1.0:
                        # 朝勤務
                        if k == 0:
                            sum_time = sum_time + 3.0
                        # 昼勤務
                        elif k == 1:
                            sum_time = sum_time + 3.0
                        # 夜勤務
                        elif k == 2:
                            sum_time = sum_time + 2.0
                        # 深夜勤務
                        elif k == 3:
                            sum_time = sum_time + 8.0
            
            # 朝昼夜深夜の時間帯で4つの個体で表す。
            # 8時間よりも多い勤務はペナルティを加える
            if sum_time > 8.0:
                penalty = penalty + 10.0
    return penalty

def select(ga, elite):
    """選択関数です。エリート選択を行います
    評価が高い順番にソートを行った後、一定以上
    :param ga: 選択を行うgenomClassの配列
    :return: 選択処理をした一定のエリート、genomClassを返す
    """
    # 現行世代個体集団の評価を低い順番にソートする
    sort_result = sorted(ga, reverse=False, key=lambda u: u.evaluation)
    # 一定の上位を抽出する
    result = [sort_result.pop(0) for i in range(elite)]
    return result

def crossover(ga_one, ga_second):
    """交叉関数です。二点交叉を行います。
    :param ga: 交叉させるgenomClassの配列
    :param ga_one:
    :param ga_second:
    :return: 二つの子孫genomClassを格納したリスト返す
    """
    # 子孫を格納するリストを生成します
    genom_list = []
    # 入れ替える二点の点を設定します
    cross_one = random.randint(0, int(GENOM_LENGTH / 4 / DAY))
    cross_second = random.randint(cross_one, int(GENOM_LENGTH / 4 / DAY))
    # 遺伝子を取り出します
    one = ga_one.getGenom()
    second = ga_second.getGenom()
    
    one_genome_list = []
    second_genome_list = []

    one_genome_all_list = []
    second_genome_all_list = []

    convert_one_genome_list = []
    convert_second_genome_list = []

    # 交叉させる
    # N人数分の遺伝子
    for i in range(len(one)):
        # 一人のDAY分のゲノムをフラット化する
        convert_one_genome_list.append(sum(one[i], []))
    for i in range(len(second)):
        # 一人のDAY分のゲノムをフラット化する
        convert_second_genome_list.append(sum(second[i], []))

    #print("conv_one",convert_one_genome_list)
    for one2, second2 in zip_longest(convert_one_genome_list, convert_second_genome_list):
        one_genome_list.append(one2[:cross_one] + second2[cross_one:cross_second] + one2[cross_second:])
        second_genome_list.append(second2[:cross_one] + one2[cross_one:cross_second] + second2[cross_second:])

    # genomClassインスタンスを生成して子孫をリストに格納する
    
    for i in range(len(one_genome_list)):
        one_genome_all_list.append(np.reshape(one_genome_list[i], (DAY, 4)).tolist())
    for i in range(len(second_genome_list)):
        second_genome_all_list.append(np.reshape(second_genome_list[i], (DAY, 4)).tolist())
    genom_list.append(Genom(one_genome_all_list, 0.0))
    genom_list.append(Genom(second_genome_all_list, 0.0))
    return genom_list


def next_generation_gene_create(ga, ga_elite, ga_progeny):
    """
    世代交代処理を行います
    :param ga: 現行世代個体集団
    :param ga_elite: 現行世代エリート集団
    :param ga_progeny: 現行世代子孫集団
    :return: 次世代個体集団
    """
    # 現行世代個体集団の評価を低い順番にソートする
    next_generation_geno = sorted(ga, reverse=False, key=lambda u: u.evaluation)
    # 追加するエリート集団と子孫集団の合計ぶんを取り除く
    [next_generation_geno.pop(0) for i in range(0, len(ga_elite) + len(ga_progeny))]
    # エリート集団と子孫集団を次世代集団を次世代へ追加します
    next_generation_geno.extend(ga_elite)
    next_generation_geno.extend(ga_progeny)
    return next_generation_geno


def mutation(ga, individual_mutation, genom_mutation):
    """突然変異関数です。
    :param ga: genomClass
    :return: 突然変異処理をしたgenomClassを返す"""
    ga_list = []
    for i in ga:
        # 個体に対して一定の確率で突然変異が起きる
        if individual_mutation > (random.randint(0, 100) / Decimal(100)):
            genom_list = []
            pre_genom_list = []
            pre_genom_list2 = []
            for i_ in i.getGenom():
                # 個体の遺伝子情報一つ一つに対して突然変異がおこる
                if genom_mutation > (random.randint(0, 100) / Decimal(100)):
                    for k in range(len(i_)):
                        for l in range(4):
                            pre_genom_list.append(float(random.randint(0, 1)))
                        pre_genom_list2.append(pre_genom_list)
                        pre_genom_list = []
                    genom_list.append(pre_genom_list2)
                    pre_genom_list2 = []
                else:
                    genom_list.append(i_)
            i.setGenom(genom_list)
            ga_list.append(i)
        else:
            ga_list.append(i)
    return ga_list


if __name__ == '__main__':
    # 一番最初の現行世代個体集団を生成します。
    current_generation_individual_group = []
    for i in range(MAX_GENOM_LIST):
        current_generation_individual_group.append(create_genom())

    for count_ in range(1, MAX_GENERATION + 1):
        # 現行世代個体集団の遺伝子を評価し、genomClassに代入します
        for i in range(MAX_GENOM_LIST):
            evaluation_result = evaluation(current_generation_individual_group[i])
            current_generation_individual_group[i].setEvaluation(evaluation_result)
        # エリート個体を選択します(評価関数で評価してソートされる)
        elite_genes = select(current_generation_individual_group, SELECT_GENOM)
        # エリート遺伝子を交叉させ、リストに格納します
        progeny_gene = []
        for i in range(0, SELECT_GENOM):
            progeny_gene.extend(crossover(elite_genes[i - 1], elite_genes[i]))
        # 次世代個体集団を現行世代、エリート集団、子孫集団から作成します
        next_generation_individual_group = next_generation_gene_create(current_generation_individual_group, elite_genes, progeny_gene)
        # 次世代個体集団全ての個体に突然変異を施します。
        # (エラーが起きる)
        next_generation_individual_group = mutation(next_generation_individual_group,INDIVIDUAL_MUTATION,GENOM_MUTATION)

        # 1世代の進化的計算終了。評価に移ります

        # 各個体適用度を配列化します。
        fits = [i.getEvaluation() for i in current_generation_individual_group]
        # 進化結果を評価します
        min_ = min(fits)
        max_ = max(fits)
        # 調和平均
        # avg_ =  pow(Decimal(sum(np.power(fits , -1)) ) / Decimal(len(fits)),-1)

        # 算術平均
        avg_ = Decimal(sum(fits)) / Decimal(len(fits)) 

        # 現行世代の進化結果を出力します
        if count_  %50 == 0  or count_ == 1:
            print("-----第{}世代の結果-----".format(count_))
            print("  Min:{}".format(min_))
            print("  Max:{}".format(max_))
            print("  Avg:{}".format(avg_))

        # 現行世代と次世代を入れ替えます
        current_generation_individual_group = next_generation_individual_group

       
    # 最終結果出力
    print("最も優れた個体は{}".format(np.reshape(elite_genes[0].getGenom(), (int(GENOM_LENGTH / DAY / 4),DAY,4)).tolist()))
```

###出力

```text

-----第1世代の結果-----
  Min:730.0
  Max:1110.0
  Avg:914.5
-----第50世代の結果-----
  Min:300.0
  Max:1040.0
  Avg:566.7
-----第100世代の結果-----
  Min:200.0
  Max:1070.0
  Avg:507.8
-----第150世代の結果-----
  Min:160.0
  Max:1080.0
  Avg:472
最も優れた個体は[[[0.0, 1.0, 1.0, 0.0], [1.0, 1.0, 0.0, 0.0], [1.0, 1.0, 0.0, 0.0], [0.0, 0.0, 0.0, 1.0], [0.0, 1.0, 0.0, 0.0], [1.0, 0.0, 0.0, 1.0], [1.0, 0.0, 1.0, 0.0]], [[0.0, 1.0, 0.0, 0.0], [1.0, 1.0, 1.0, 0.0], [0.0, 0.0, 0.0, 1.0], [0.0, 0.0, 1.0, 0.0], [0.0, 1.0, 0.0, 0.0], [1.0, 0.0, 1.0, 1.0], [1.0, 0.0, 0.0, 0.0]], [[0.0, 0.0, 0.0, 1.0], [1.0, 0.0, 1.0, 0.0], [1.0, 0.0, 0.0, 0.0], [1.0, 1.0, 0.0, 0.0], [1.0, 0.0, 1.0, 1.0], [1.0, 0.0, 0.0, 0.0], [1.0, 1.0, 0.0, 0.0]], [[0.0, 0.0, 0.0, 1.0], [0.0, 0.0, 0.0, 1.0], [1.0, 1.0, 1.0, 0.0], [1.0, 0.0, 0.0, 0.0], [0.0, 1.0, 1.0, 0.0], [1.0, 1.0, 1.0, 0.0], [1.0, 0.0, 0.0, 0.0]], [[0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0], [1.0, 1.0, 1.0, 0.0], [1.0, 0.0, 0.0, 0.0], [1.0, 0.0, 1.0, 0.0], [1.0, 1.0, 1.0, 1.0], [1.0, 1.0, 1.0, 0.0]], [[0.0, 0.0, 0.0, 0.0], [1.0, 1.0, 1.0, 0.0], [1.0, 0.0, 1.0, 0.0], [0.0, 0.0, 1.0, 0.0], [0.0, 0.0, 0.0, 0.0], [0.0, 1.0, 1.0, 0.0], [0.0, 0.0, 0.0, 1.0]], [[0.0, 0.0, 1.0, 0.0], [1.0, 0.0, 1.0, 0.0], [1.0, 1.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0], [0.0, 0.0, 0.0, 1.0], [1.0, 0.0, 1.0, 1.0], [1.0, 1.0, 0.0, 0.0]], [[0.0, 0.0, 0.0, 1.0], [0.0, 1.0, 1.0, 0.0], [0.0, 1.0, 0.0, 0.0], [0.0, 1.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0], [1.0, 0.0, 0.0, 0.0], [1.0, 1.0, 1.0, 0.0]], [[1.0, 0.0, 0.0, 1.0], [0.0, 1.0, 0.0, 0.0], [1.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 1.0], [1.0, 1.0, 1.0, 0.0], [1.0, 0.0, 0.0, 0.0], [1.0, 0.0, 1.0, 0.0]], [[0.0, 0.0, 0.0, 0.0], [1.0, 0.0, 1.0, 0.0], [1.0, 0.0, 1.0, 0.0], [0.0, 0.0, 1.0, 0.0], [1.0, 0.0, 0.0, 0.0], [1.0, 1.0, 1.0, 0.0], [1.0, 1.0, 1.0, 1.0]], [[0.0, 0.0, 0.0, 0.0], [1.0, 0.0, 1.0, 0.0], [0.0, 0.0, 0.0, 0.0], [0.0, 1.0, 1.0, 0.0], [0.0, 1.0, 1.0, 0.0], [0.0, 0.0, 0.0, 1.0], [1.0, 0.0, 1.0, 0.0]], [[1.0, 1.0, 1.0, 0.0], [1.0, 1.0, 0.0, 0.0], [1.0, 1.0, 0.0, 0.0], [1.0, 0.0, 1.0, 0.0], [0.0, 1.0, 0.0, 1.0], [0.0, 0.0, 1.0, 0.0], [0.0, 0.0, 0.0, 1.0]], [[1.0, 0.0, 1.0, 0.0], [1.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 1.0], [0.0, 0.0, 0.0, 1.0], [0.0, 0.0, 1.0, 0.0], [0.0, 0.0, 0.0, 1.0], [1.0, 0.0, 1.0, 0.0]], [[1.0, 1.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0], [1.0, 0.0, 0.0, 0.0], [1.0, 1.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0], [0.0, 1.0, 0.0, 0.0], [0.0, 1.0, 0.0, 0.0]], [[0.0, 0.0, 0.0, 1.0], [0.0, 0.0, 0.0, 1.0], [1.0, 1.0, 0.0, 0.0], [0.0, 1.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0], [0.0, 0.0, 1.0, 0.0], [1.0, 1.0, 1.0, 0.0]], [[1.0, 1.0, 0.0, 0.0], [1.0, 1.0, 1.0, 0.0], [0.0, 1.0, 0.0, 0.0], [0.0, 1.0, 1.0, 0.0], [0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0], [0.0, 0.0, 0.0, 1.0]], [[1.0, 1.0, 1.0, 0.0], [0.0, 1.0, 0.0, 0.0], [0.0, 1.0, 0.0, 1.0], [0.0, 1.0, 1.0, 0.0], [0.0, 1.0, 1.0, 0.0], [0.0, 1.0, 1.0, 0.0], [0.0, 0.0, 0.0, 0.0]], [[0.0, 0.0, 0.0, 1.0], [1.0, 0.0, 0.0, 0.0], [0.0, 1.0, 0.0, 0.0], [0.0, 1.0, 1.0, 0.0], [1.0, 1.0, 1.0, 0.0], [0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 1.0, 1.0]], [[1.0, 0.0, 1.0, 0.0], [0.0, 1.0, 1.0, 0.0], [0.0, 1.0, 1.0, 0.0], [0.0, 1.0, 0.0, 1.0], [1.0, 1.0, 1.0, 0.0], [0.0, 0.0, 0.0, 0.0], [1.0, 1.0, 1.0, 0.0]], [[1.0, 1.0, 0.0, 0.0], [0.0, 0.0, 0.0, 1.0], [0.0, 1.0, 0.0, 0.0], [1.0, 1.0, 0.0, 0.0], [1.0, 1.0, 0.0, 0.0], [0.0, 1.0, 1.0, 0.0], [0.0, 0.0, 1.0, 0.0]], [[0.0, 0.0, 1.0, 1.0], [0.0, 0.0, 1.0, 0.0], [0.0, 0.0, 0.0, 1.0], [1.0, 1.0, 0.0, 0.0], [0.0, 0.0, 0.0, 1.0], [0.0, 1.0, 0.0, 0.0], [1.0, 0.0, 1.0, 0.0]], [[0.0, 1.0, 0.0, 0.0], [1.0, 1.0, 0.0, 1.0], [1.0, 0.0, 0.0, 0.0], [1.0, 0.0, 0.0, 0.0], [1.0, 1.0, 0.0, 1.0], [0.0, 0.0, 0.0, 0.0], [0.0, 1.0, 1.0, 0.0]], [[1.0, 1.0, 0.0, 1.0], [1.0, 1.0, 0.0, 0.0], [1.0, 1.0, 1.0, 0.0], [0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 1.0], [0.0, 1.0, 0.0, 0.0], [0.0, 1.0, 1.0, 0.0]], [[0.0, 0.0, 1.0, 0.0], [0.0, 1.0, 0.0, 1.0], [0.0, 0.0, 1.0, 0.0], [0.0, 1.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0], [1.0, 1.0, 0.0, 0.0], [1.0, 0.0, 1.0, 0.0]], [[1.0, 1.0, 0.0, 0.0], [1.0, 0.0, 0.0, 0.0], [1.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0], [0.0, 1.0, 1.0, 0.0], [1.0, 1.0, 1.0, 0.0], [0.0, 1.0, 0.0, 0.0]], [[0.0, 0.0, 0.0, 0.0], [0.0, 1.0, 0.0, 0.0], [1.0, 0.0, 0.0, 0.0], [1.0, 1.0, 0.0, 0.0], [0.0, 0.0, 0.0, 1.0], [0.0, 1.0, 0.0, 0.0], [1.0, 1.0, 0.0, 0.0]], [[0.0, 0.0, 0.0, 1.0], [0.0, 0.0, 0.0, 0.0], [1.0, 1.0, 1.0, 0.0], [0.0, 0.0, 0.0, 1.0], [0.0, 1.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0], [1.0, 0.0, 0.0, 0.0]], [[0.0, 0.0, 0.0, 1.0], [0.0, 0.0, 0.0, 0.0], [1.0, 0.0, 0.0, 0.0], [0.0, 1.0, 0.0, 0.0], [1.0, 1.0, 0.0, 0.0], [1.0, 1.0, 0.0, 0.0], [0.0, 0.0, 0.0, 1.0]], [[1.0, 1.0, 0.0, 0.0], [1.0, 1.0, 1.0, 0.0], [1.0, 0.0, 1.0, 0.0], [0.0, 0.0, 0.0, 1.0], [0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 1.0]], [[1.0, 0.0, 1.0, 0.0], [0.0, 1.0, 1.0, 0.0], [0.0, 1.0, 1.0, 0.0], [1.0, 0.0, 1.0, 0.0], [0.0, 0.0, 0.0, 1.0], [0.0, 1.0, 1.0, 0.0], [0.0, 1.0, 1.0, 0.0]]]
```