第7回: アニーリング基礎6.補助スピン活用基礎2（シフト最適化問題）

-5人の従業員の14日分のシフトの組み合わせを、生産性と制約条件を加味しながら算出する。
前回：第6回: アニーリング基礎5.補助スピン活用基礎1（ナップザック問題）とTTSの活用
Keyword：*シフト最適化問題、補助スピン、TTS

1. 基本準備/設定

1.1 量子アニーリング用モジュールのインストール

まず、初めに、以下のモジュールをインストールする。

• dwave-ocean-sdk：Dwaveアニーリングマシンのモジュール
• openjij：アニーリングシュミレーター
• pyqubo：Qubo行列作成サポートモジュール
• optuna：ハイパーパラメタの最適化

pip install dwave-ocean-sdk
pip install openjij
pip install pyqubo
pip install optuna

1.2 Dwaveの初期設定

続いて、Dwaveのtoken、enpoint(アクセスURL)、使用するマシンの指定まで実施。併せて、Jijシュミレーターも設定する。

###基本設定
token = "DEV-73ac79954ef21XXXXXXXXXXXXXXXXXXXXXX"
endpoint = "https://cloud.dwavesys.com/sapi/"
solver = "Advantage2_prototype1.1"
###Dwaveの設定
from dwave.system import DWaveSampler, EmbeddingComposite
dw_sampler = DWaveSampler(solver=solver, token=token, endpoint = endpoint)
###Jijシュミレーターの設定
from openjij import SQASampler
sampler_jij = SQASampler()

1.3 基本ライブラリのインポート

import numpy as np
import pandas as pd

2. 問題設定

今回、$N=5$人の従業員の、$D=14$日分のシフトを作成する。

work_member = ["Aさん","Bさん","Cさん","Dさん","Eさん"]
work_day = ["月","火","水","木","金","土","日","月","火","水","木","金","土","日"]
M = len(work_member)
D = len(work_day)

シフトは以下の基本ルール(=制約条件)を満たす。

• 制約1: 平日の最低出勤数は3人
• 制約2: 休日の最低出勤数は2人
• 制約3: 最低勤務日数は8日
• 制約4: 全員1回以上は土日に休む
• 制約5: 3連休は禁止

また、追加情報として、表1.シフト希望、表2.従業員の生産性、表3.従業員同士の相性が与えられる。
それぞれの従業員の生産性は、それぞれの生産性$×$その他全従業員との相性合計により計算する。
本問題では、制約条件を満たしつつ、シフト希望に沿いつつ、14日間全体の生産性を最大化させることを最適化の目的とする。
表1.シフト希望は以下とする。

shift_A, shift_B, shift_C, shift_D, shift_E = ["〇","〇","〇","〇","×","〇","△","×","×","〇","△","〇","〇","×"], ["×","×","×","〇","〇","〇","△","〇","×","〇","〇","〇","〇","〇"] ,["△","△","△","×","〇","×","〇","×","〇","〇","〇","△","△","△"] ,["〇","〇","×","×","×","×","〇","〇","〇","〇","△","×","×","×"] ,["〇","〇","△","×","×","〇","×","〇","〇","×","×","〇","△","△"]
shift_original = np.array([shift_A,shift_B,shift_C,shift_D,shift_E])
df_shift_original=pd.DataFrame(data = shift_original,columns = work_day, index = work_member)
df_shift_original

後続処理の為、「〇」=0、「△」=1、「×」=2として、数値に変換する。

S = np.where(shift_original=="〇",0,np.where(shift_original=="△",1,2))#Shift
#df_shift=pd.DataFrame(data = S,columns = work_day, index = work_member)
#df_shift

表2.従業員の生産性は以下とする。

P = [10,7,5,3,2] #productivity
df_prod=pd.DataFrame(data=P,index=work_member)
df_prod

表3.従業員の相性は以下とする。

C = np.array([[0,-2,5,2,1],[-2,0,-3,2,5],[3,-3,0,-3,-4],[2,2,-3,0,5],[1,5,-4,5,0]])#Compatibility
df_comp=pd.DataFrame(data = C,columns=work_member,index=work_member)
df_comp #Compatibility

3. 定式化とQUBO行列の取出し

3.0 スピンの設定

今回、設定するスピンは$x_{md}$として、例えば$x_{23}$は$m=2$のCさんが、$d=3$の木に働く場合は1、働かない場合を0として設定する。尚、今回の定式化に於いては全て項に$\lambda$を設定する。

from pyqubo import Array, Placeholder
x = Array.create(name="x", shape =(M,D), vartype="BINARY") 

3.1 効率性Pの最大化

効率性Pは、従業員の効率性×その従業員と他の出勤従業員との相性によって決定される。効率性は最大化対象なので、符号はマイナスとなる。つまり、以下のように表現できる。

E_{1}(x) = －\lambda_{e1}\sum_{d=1}^D \sum_{m=1}^M \sum_{n=1}^M P_{m} x_{md} x_{nd} C_{mn}

E1 = 0
for d in range(D):
  for m in range(M):
    for n in range(M):
      E1 = E1 - P[m]*x[m,d]*x[n,d]*C[m,n] #λは後で設定する

3.2 シフト希望Sの充足

シフト希望$S$は「〇」=0、「△」=1、「×」=2で定義した。仕事をする場合、$x_{md}=1$となるが、「〇」の場合に業務行えば不満にならない。一方で「×」の場合に働けば不満につながる。下式のように$x_{md}$とシフト希望日を積すれば不満度を表現できる。尚、業務しない場合は、$x_{md}=0$となり、積しても0のままであるが、働かないので不満は発生しないとの考え方である。尚、これは不満度なので、コスト関数である。

C_{1}(x) = \lambda_{c1}\sum_{m=1}^M \sum_{d=1}^D x_{md}S_{md}

C1=0
for m in range(M):
  for d in range(D):
    C1 = C1 + x[m,d]*S[m,d] #λは後で設定する

3.3 制約1: 平日の最低出勤数は3人

出勤者数は3≦出勤者数≦5となり、独立型の補助スピン$y_{k}$を用いる場面であるが、スピンが増大することを鑑みて、$x_{md}$のみで対応する。3～5人が出勤すればよい為、その中間をとって4人が業務するようにコスト関数を追加する。

C_{2}(x) = \lambda_{c2}\sum_{d=1}^{D_{wd}} (\sum_{m=1}^M x_{md}-4)^2\\
ここで、D_{wd}は平日を意味する。

C2=0
for d in range(D):
  if not d == 5 and not d ==6 and not d ==12 and not d ==13:
    _temp = 0
    for m in range(M):
      _temp = _temp +(x[m,d]-4)**2
  C2 = C2 + _temp

3.4 制約2: 休日の最低出勤数は2人

考え方は上の制約1と同じ。2～5人が出勤すればよい為、その中間をとって3.5人が業務するようにコスト関数を追加する。

C_{3}(x) = \lambda_{c3}\sum_{d=1}^{D_{we}} (\sum_{m=1}^M x_{md}-3.5)^2\\
ここで、D_{we}は週末を意味する。

C3=0
for d in range(D):
  if d == 5 and d ==6 and d ==12 and d ==13:
    _temp = 0
    for m in range(M):
      _temp = _temp +(x[m,d]-3.5)**2
  C3 = C3 + _temp

3.5 制約3: 最低勤務日数は8日

これも上の制約と同じ。8～14の中間の11日業務するように設定すればよい。

C_{4}(x) = \lambda_{c4}\sum_{m=1}^{M} (\sum_{d=1}^D x_{md}-11)^2

C4=0
for m in range(M):
  _temp = 0
  for d in range(D):
    _temp = _temp +(x[m,d]-11)**2
  C4 = C4 + _temp

3.6 制約4: 全員1回以上は土日に休む

1～4にするという不等式なので本来は補助スピンを使用すべきだが、Qubo行列を小さくするために、これを上の制約条件と同じように実施する。土日は全てで4回あるので、1～4の中間の2.5と設定できる。しかし、その場合、土日休みが多くなりすぎるため、1～4の内、できる限り1に近いとして2.0とする。

C_{5}(x) = \lambda_{c5}\sum_{m=1}^{M} (\sum_{d=1}^{D_{we}} x_{md}-2)^2

C5 = 0
for m in range(M):
  for d in range(D):
    if d == 5 and d ==6 and d ==12 and d ==13:
      _temp = 0
      _temp = _temp + (x[m,d]-2)**2
  C5 = C5 + _temp

3.7 制約5: 3連休は禁止

3回連続で休むとコストが増加するように以下のコスト関数を使用する。

C_{6}(x) = \lambda_{c6}\sum_{m=1}^{M} \sum_{d=1}^{D-2} (1-x_{md})(1-x_{md+1})(1-x_{md+2})

この式は3次元なのでアニーリングでは処理できない。そこで、2次元にするべく補助スピンを導入する。今回は$x_{md}$と$y_{k}$が関連する関連型。この時、どういう関係式にするか定義が必要。そこで上の式を展開する。

C_{6}(x) = \lambda_{c6}\sum_{m=1}^{M} \sum_{d=1}^{D-2} (1-x_{md}-x_{md+1}+x_{md}x_{md+1})(1-x_{md+2})

ここで、$x_{md}x_{md+1}$に$x_{md+2}$を積することで、この部分だけが3次式となる。
そこで、$y_{k}$を以下とする。

y_{md} = x_{md}x_{md+1}

元の式に代入して、

C_{6a}(x) = \lambda_{c6}\sum_{m=1}^{M} \sum_{d=1}^{D-2} (1-x_{md}-x_{md+1}+y_{md})(1-x_{md+2})

ここで、$y_{md} = x_{md}x_{md+1}$の関係を担保するために、以下の制約式を追加する必要がある。この制約式の形になる理由は後述。

C_{6b}(x) = \lambda_{c6}\sum_{m=1}^{M} \sum_{d=1}^{D-2}(3y_{md}+x_{md}x_{md+1}-2y_{md}x_{md}-2y_{md}x_{md+1})

尚、$y_{md}$を見ると、mは1～Mまで、dは1～D-2までに範囲を取るため、$M×(D-2)$の行が増大することになる。今回$M=5、D=14$なので、$x_{md}$の個数は$5×14=70$、$y_{md}$の個数は、$5×12=60$となるので、今回のQubo行列は$70×70=4900$→$130×130=16900$まで増加することになる。

### 補助スピンの設定
from pyqubo import Array, Placeholder
y = Array.create(name="y", shape =(M,D-2), vartype="BINARY") 

###コスト関数C5の定式化
C6a =0
for m in range(M):
  for d in range(D-2):
    C6a = C6a +(1-x[m,d]-x[m,d+1]+y[m,d])*(1-x[m,d+2])

###コスト関数C5の補助スピンを担保する制約式の定式化
C6b = 0
for m in range(M):
  for d in range(D-2):
    C6b = C6b + (3*y[m,d]+x[m,d]*x[m,d+1]-2*y[m,d]*x[m,d]-2*y[m,d]*x[m,d+1])

3.8 ハミルトニアンの作成

上記で整理してきたエネルギー式、コスト関数をまとめる。

H = \lambda_{e1}E_{1}(x)+
    \lambda_{c1}C_{1}(x)+
    \lambda_{c2}C_{2}(x)+
    \lambda_{c3}C_{3}(x)+
    \lambda_{c4}C_{4}(x)+
    \lambda_{c5}C_{5}(x)+
    \lambda_{c6}(C_{6a}(x)+C_{6b}(x))

### λの設定
from pyqubo import Placeholder
lambda_e1,lambda_c1,lambda_c2,lambda_c3,lambda_c4,lambda_c5,lambda_c6 = Placeholder("lambda_e1"),Placeholder("lambda_c1"),Placeholder("lambda_c2"), Placeholder("lambda_c3"),Placeholder("lambda_c4"),Placeholder("lambda_c5"),Placeholder("lambda_c6")

### ハミルトニアンの作成
H = lambda_e1*E1 + lambda_c1*C1 + lambda_c2*C2 + lambda_c3*C3 + lambda_c4*C4 + lambda_c5*C5 + lambda_c6*(C6a + C6b)

4. 計算の実行

4.1 λの設定 下準備:各ハミルトニアンの可視化

まず、ハミルトニアンを構成する各ハミルトニアンの大きさ、特徴を確認する。さらに、制約条件の優先度を鑑みて、λの範囲をある程度絞る。まずPyquboの辞書型Quboをリスト型に変換してヒートマップで可視化する関数を定義する。

import seaborn as sns
import matplotlib.pyplot as plt

def transform_qubo(_H,_m,_d, _show = True,_size=5):

  ###QUBO行列の取出し(ただし辞書型)、行列に置き換える
  _qubo_box = np.zeros(((_m*_d)+(_m*(_d-2)))**2).reshape((_m*_d)+(_m*(_d-2)),(_m*_d)+(_m*(_d-2))) 
  _qubo = _H.compile().to_qubo(feed_dict={})[0]

  for _key,_value in zip(list(_qubo.keys()),list(_qubo.values())):
    
    ###x(列)が何列目か計算
    _temp_x = _key[0].replace("]","").split("[") #"x[3][4]"から]をreplaceで削除して"x[3[4"に変換。続いて、[でsplitして["x","3","4"]のリストに変換。
    #print(_temp_x)
    if _temp_x[0] =="x":
      _x = int(_temp_x[1])*D+int(_temp_x[2])
    elif _temp_x[0] =="y":
      _x = int(_temp_x[1])*(D-2)+int(_temp_x[2])+M*D
    
    ###y(行)が何列目か計算
    _temp_y = _key[1].replace("]","").split("[") #"x[3][4]"から]をreplaceで削除して"x[3[4"に変換。続いて、[でsplitして["x","3","4"]のリストに変換。
    #print(_temp_y)
    if _temp_x[0] =="x":
      _y = int(_temp_x[1])*D+int(_temp_x[2])
    elif _temp_x[0] =="y":
      _y = int(_temp_x[1])*(D-2)+int(_temp_x[2])+M*D
    
    ###行列に変換する
    _qubo_box[_x][_y] = _value

    #可視化
  if _show == True:
    _qubo_mask = (_qubo_box==0)
    plt.figure(figsize = (_size,_size))
    return sns.heatmap(_qubo_box, square=True,robust=False,linewidths=0.0005,linecolor="aliceblue",cmap="Reds",mask = _qubo_mask)  

  else:
    return _qubo_box

print("E1:従業員の生産性"),transform_qubo(E1,M,D)

print("C1: シフト希望S充足度"),transform_qubo(C1,M,D)

print("C2: 制約1: 平日の最低出勤数は3人"),transform_qubo(C2,M,D)

print("C3: 制約2: 休日の最低出勤数は2人"),transform_qubo(C3,M,D)

print("C4: 制約3: 最低勤務日数は8日"),transform_qubo(C4,M,D)

print("C5: 制約4: 全員1回以上は土日に休む"),transform_qubo(C5,M,D)

print("C6: 制約5: 3連休は禁止"),transform_qubo(C6a,M,D),transform_qubo(C6b,M,D)

結局λを決定するにあたって、それぞれの最大値がおおむね同じぐらいの値を示せばよい。そこで、それぞれのハミルトニアンの最大値を取得して、その最大にするための倍率を取得する。これを用いてハミルトニアンの定式化を実施する。

max_list = [np.amax(np.abs(transform_qubo(_target,M,D,False))) for _target in [E1,C1,C2,C3,C4,C5,C6a,C6b]]
lambda_guide = [round(1/_temp,2) for _temp in list(max_list / max(max_list))]
lambda_guide

<Output>
[1.62, 52.5, 5.0, 1.07, 5.0, 1.0, 35.0, 35.0]

4.2 定式化

先ほど算出のλガイドを参考にQubo行列を抜き出し、Dwaveで計算を行う。

Qubo = H.compile().to_qubo\
        (feed_dict=\
         {"lambda_e1":2,\
          "lambda_c1":200,\
          "lambda_c2":5,\
          "lambda_c3":1.5,\
          "lambda_c4":5,\
          "lambda_c5":1,\
          "lambda_c6":100})[0]
sampler = EmbeddingComposite(dw_sampler)
solution = sampler.sample_qubo(Qubo, num_reads=100)

結果を可視化する関数を定義する。

def evaluation(_solution,_m,_d,_e):
  #####データ変換
  ###Dwave出力生データをシフト表形式に変換
  _result = np.zeros(_m*_d).reshape(_m,_d)
  for _value,_target in zip(solution.record[_e][0],solution.variables):
    _target = _target.replace("]","").split("[")
    if _target[0][0] == "x":
      _result[int(_target[1])][int(_target[2])] = _value
  _result

  #####Evaluation
  ###制約1、2: 平日週末の最低出勤数の確認
  ##3人以上かどうか確認
  if sum(np.where(sum(_result)>=3,0,1)) == 0:
    j12 = 1
  ##週末、2人以上かどうか確認
  else:
    if np.dot(np.where(sum(_result)>=2,0,1),[0,0,0,0,0,1,1,0,0,0,0,0,1,1]) == 0:
      j12 = 1
    else:
      j12 = 0

  ###制約3: 各従業員の出勤数を確認
  if sum(np.where(np.sum(_result,axis=1)>=8,0,1)) == 0:
    j3 = 1
  else:
    j3 = 0

  ###制約4: 全員1回以上は土日に休む
  we = [0,0,0,0,0,1,1,0,0,0,0,0,1,1]
  if sum(np.where(np.array([np.dot(__temp,we) for __temp in _result])>=1,0,1)) == 0:
    j4 = 1
  else:
    j4 = 0

  ###制約5: 3連休禁止
  _judge = [0 for abc in range(M)]
  for _k,_temp in enumerate(_result):
    _test = 0
    for __temp in _temp:
      if __temp == 0:
        _test = _test + 1
        if _test == 3:
          _judge[_k] = 1
      else:
        _test = 0
  if sum(_judge) == 0:
    j5 = 1
  else:
    j5 = 0

  ###効率性の計算
  _E = 0
  for d in range(_d):
    for m in range(_m):
      for n in range(_m):
        _E = _E +P[m]*_result[m,d]*_result[n,d]*C[m,n] 

  ###シフト希望充足率
  _C=0
  for m in range(_m):
    for d in range(_d):
      _C = _C + _result[m,d]*S[m,d] #λは後で設定する

  ### Output
  return [j12,j3,j4,j5,_E,round(1/_C,2),_result]

#evaluation(solution,M,D,target)

Evaluation関数を用いて結果を可視化する。出力は左から順に[制約式12,制約式3,制約式4,制約式5,生産性,シフト充足率]を意味する。試しに10件の計算結果を確認してみる。

for __k in range(10):
  print(evaluation(solution,M,D,__k)[:-1])

<Output>
[1, 0, 1, 0, 692.0, 0.06]
[1, 1, 1, 0, 532.0, 0.07]
[1, 0, 1, 0, 566.0, 0.06]
[1, 1, 1, 0, 491.0, 0.08]
[1, 1, 1, 0, 701.0, 0.07]
[1, 1, 1, 0, 753.0, 0.06]
[1, 1, 1, 0, 749.0, 0.05]
[1, 0, 1, 1, 526.0, 0.06]
[1, 1, 1, 0, 659.0, 0.06]
[1, 1, 1, 0, 567.0, 0.07]

一番良い結果を抽出してみる。

_result_box = []
###Evaluation関数の全ての結果を積してリスト化する
for __k in range(len(solution.record)):
  _result = 1
  for _temp in evaluation(solution,M,D,__k)[:-1]:
    _result = _result*_temp
  _result_box.append(_result)
###リストの中の最大値のIndexを呼び出す
target = _result_box.index(max(_result_box))
#print(f"n={target}, {evaluation(solution,M,D,target)[:-1]}")

df_result = pd.DataFrame(data = np.where(evaluation(solution,M,D,target)[-1]==1,"〇","×"),index = work_member, columns = work_day)
print(evaluation(solution,M,D,target)[:-1])
df_result

df_shift_original

概ね制約条件を満たしつつ、シフト希望を満たしていることが確認できた。続いて、ハイパーパラメタの最適化を行う。

4.3 ハイパーパラメタの最適化

設定したλの0.75倍～5倍の範囲に設定して最適なλの組み合わせを探索してみる。

###ライブラリの呼び出し
import optuna

###アニーリング計算実行関数
def a_cal(_Qubo, machine = "Dwave", case =100):
### Dwaveかどうかの選択
  if machine == "Dwave":
    sampler_dw = EmbeddingComposite(dw_sampler)
    _solution = sampler_dw.sample_qubo(_Qubo, num_reads=case)
    
### Jijかどうかの選択
  elif machine =="Jij":
    sampler_jij = SQASampler()
    _solution = sampler_jij.sample_qubo(_Qubo,num_reads=case)
  return(_solution)

###関数の定義
def objective(trial):

  ###ハイパーパラメタの設定
  _lambda_e1, _lambda_c1, _lambda_c2, _lambda_c3, _lambda_c4, _lambda_c5, _lambda_c6,=\
        trial.suggest_uniform('_lambda_e1', 1.5, 2.5),\
        trial.suggest_uniform('_lambda_c1', 150, 250),\
        trial.suggest_uniform('_lambda_c2', 3.75, 6.25),\
        trial.suggest_uniform('_lambda_c3', 1.125, 1.875),\
        trial.suggest_uniform('_lambda_c4', 3.75, 6.25),\
        trial.suggest_uniform('_lambda_c5', 0.75, 1.25),\
        trial.suggest_uniform('_lambda_c6', 75, 125)

  ###コンパイルとQuboの取り出し
  _model = H.compile()
  _Qubo, offset = _model.to_qubo(feed_dict=\
                                 {"lambda_e1":_lambda_e1,\
                                  "lambda_c1":_lambda_c1,\
                                  "lambda_c2":_lambda_c2,\
                                  "lambda_c3":_lambda_c3,\
                                  "lambda_c4":_lambda_c4,\
                                  "lambda_c5":_lambda_c5,\
                                  "lambda_c6":_lambda_c6})

  ###a_calでアニーリング計算実施し、それをrateで正答率にする。
  _scoure = 0
  _solution = a_cal(_Qubo, "Jij")
  for _i in range(len(_solution.record)):
    _temp = evaluation(_solution,M,D,_i)[:-1]
    #print(_temp)
    _scoure = _scoure + _temp[0]*_temp[1]*_temp[2]*_temp[3]*_temp[4]*_temp[5]
    _scoure = -1*_scoure
  return _scoure

### Optunaでの最適化
study = optuna.create_study()
study.optimize(objective, n_trials=10)

### Optunaでの最適解を表示
best_lambda_e1, best_lambda_c1, best_lambda_c2, best_lambda_c3, best_lambda_c4, best_lambda_c5, best_lambda_c6 =\
   [i for i in study.best_params.values()][0],\
   [i for i in study.best_params.values()][1],\
   [i for i in study.best_params.values()][2],\
   [i for i in study.best_params.values()][3],\
   [i for i in study.best_params.values()][4],\
   [i for i in study.best_params.values()][5],\
   [i for i in study.best_params.values()][6]
   
print(f"l_e1: {best_lambda_e1}, l_c1: {best_lambda_c1}, l_c2: {best_lambda_c2}, l_c3: {best_lambda_c3}, l_c4: {best_lambda_c4}, l_c5: {best_lambda_c5}, l_c6: {best_lambda_c6}")

<Output>
[I 2022-12-04 12:19:04,349] A new study created in memory with name: no-name-b4518aef-<ipython-input-99-77bd439ec5a3>:22: FutureWarning: suggest_uniform has been deprecated in v3.0.0. This feature will be :func:`~optuna.trial.Trial.suggest_float` instead.
  trial.suggest_uniform('_lambda_e1', 1.5, 2.5),\
[I 2022-12-04 12:19:09,745] Trial 1 finished with value: 87.99000000000001 and parameters: {'_lambda_e1': 1.8448773734221344, '_lambda_c1': 180.81584517467016, '_lambda_c2': 5.929699389401177, '_lambda_c3': 1.2795362300710176, '_lambda_c4': 4.04920168881826, '_lambda_c5': 1.0397748951204826, '_lambda_c6': 94.33718743006153}. Best is trial 0 with value: 87.99000000000001.
[I 2022-12-04 12:19:12,363] Trial 2 finished with value: 87.99000000000001 and parameters: {'_lambda_e1': 1.8981552576387324, '_lambda_c1': 190.32862746800987, '_lambda_c2': 5.89500081190968, '_lambda_c3': 1.5602585507397548, '_lambda_c4': 5.274906176132985, '_lambda_c5': 0.936517010084013, '_lambda_c6': 86.66713120484098}. Best is trial 0 with value: 87.99000000000001.
[I 2022-12-04 12:19:15,021] Trial 3 finished with value: 87.99000000000001 and parameters: {'_lambda_e1': 1.9275895427752383, '_lambda_c1': 206.88298574168613, '_lambda_c2': 4.969764459110104, '_lambda_c3': 1.8728744203408163, '_lambda_c4': 4.704927978923892, '_lambda_c5': 0.7749598383014773, '_lambda_c6': 96.66907032135344}. Best is trial 0 with value: 87.99000000000001.
[I 2022-12-04 12:19:30,801] Trial 9 finished with value: 87.99000000000001 and parameters: {'_lambda_e1': 1.8412015692070915, '_lambda_c1': 245.84917858556722, '_lambda_c2': 4.7585901485080235, '_lambda_c3': 1.7692533568594393, '_lambda_c4': 3.836999384632095, '_lambda_c5': 0.8600410677024013, '_lambda_c6': 88.62908051519719}. Best is trial 0 with value: 87.99000000000001.
l_e1: 1.7358691909513908, l_c1: 196.06570259981928, l_c2: 4.094233799761867, l_c3: 1.6960408697850107, l_c4: 4.314505465237894, l_c5: 0.8133512003305898, l_c6: 76.28482239523484

得られた最適化なパラメタを用いて結果を確認する。

Qubo = H.compile().to_qubo\
        (feed_dict=\
         {"lambda_e1":best_lambda_e1,\
          "lambda_c1":best_lambda_c1,\
          "lambda_c2":best_lambda_c2,\
          "lambda_c3":best_lambda_c3,\
          "lambda_c4":best_lambda_c4,\
          "lambda_c5":best_lambda_c5,\
          "lambda_c6":best_lambda_c6})[0]
sampler = EmbeddingComposite(dw_sampler)
solution = sampler.sample_qubo(Qubo, num_reads=100)

_result_box = []
###Evaluation関数の全ての結果を積してリスト化する
for __k in range(len(solution.record)):
  _result = 1
  for _temp in evaluation(solution,M,D,__k)[:-1]:
    _result = _result*_temp
  _result_box.append(_result)
###リストの中の最大値のIndexを呼び出す
target = _result_box.index(max(_result_box))
print(f"n={target}, {evaluation(solution,M,D,target)[:-1]}")

<Output>
n=89, [1, 1, 1, 1, 665.0, 0.07]

df_result = pd.DataFrame(data = np.where(evaluation(solution,M,D,target)[-1]==1,"〇","×"),index = work_member, columns = work_day)
print(evaluation(solution,M,D,target)[:-1])
df_result

シフト希望と比較してみると以下の通り。

df_shift_original

ある程度希望にそってシフトが作成できたことを確認した。

5. 関連型の補助スピンの制約式の根拠

3回連休不可のコスト関数を設定した時、以下のようにと整理した。

C_{6}(x) = \lambda_{c6}\sum_{m=1}^{M} \sum_{d=1}^{D-2} (1-x_{md})(1-x_{md+1})(1-x_{md+2})\\
C_{6}(x) = \lambda_{c6}\sum_{m=1}^{M} \sum_{d=1}^{D-2} (1-x_{md}-x_{md+1}+x_{md}x_{md+1})(1-x_{md+2})\\
ここで、y_{md} = x_{md}x_{md+1}\\
C_{6a}(x) = \lambda_{c6}\sum_{m=1}^{M} \sum_{d=1}^{D-2} (1-x_{md}-x_{md+1}+y_{md})(1-x_{md+2})\\
C_{6b}(x) = \lambda_{c6}\sum_{m=1}^{M} \sum_{d=1}^{D-2}(3y_{md}+x_{md}x_{md+1}-2y_{md}x_{md}-2y_{md}x_{md+1})

ここでは、関連型の補助スピンを使用する時、追加の制約式が上記$C_{6b}$となる理由を整理する。

つまり、$x,y,z$のスピンがあるとして、それぞれの関係性は以下で表現できる。

R = r_{0} + r_{1}x + r_{2}y + r_{3}z + r_{4}xy + r_{5}xz + r_{6}yz + r_{7}xyz\\
ここで、z = xy　として、\\
R = r_{0} + r_{1}x + r_{2}y + r_{3}z + r_{4}xy + r_{5}xz + r_{6}yz + r_{7}z^2\\
z^2は関係性を示すものでないので除外\\
R = r_{0} + r_{1}x + r_{2}y + r_{3}z + r_{4}xy + r_{5}xz + r_{6}yz

ここで、$z=xy$で、それぞれ0or1のスピンなので、$x,y,z$がとりうるパターンは下表左、一方で右はとりえないパターンである。

$x$	$y$	$z$	$R$		$x_{error}$	$y_{error}$	$z_{error}$	$R$
0	0	0	0		0	0	1	A
0	1	0	0		0	1	1	B
1	0	0	0		1	0	1	C
1	1	1	0		1	1	0	D

とりえないパターンの時の$R$の数値をそれぞれ$A, B, C, D$として、それ以外のパターンの正しい時の$R$を$0$として、$A, B, C, D$を計算すると$A=3×n, B=C=D=n$となる。ここで$n=1$と固定すると、$r_{0}=r_{1}=r_{2}=0, r_{3}=3, r_{4}=1, r_{5}=-2, r_{6}=-2$となる。つまり、

R = r_{0} + r_{1}x + r_{2}y + r_{3}z + r_{4}xy + r_{5}xz + r_{6}yz\\
R = 0 + 0x + 0y  +3z  +1xy -2xz -2yz

このRが0でないならば、$z=xy$が成立していないことを意味する。アニーリングでは最小化の過程でこれが0になり、0にならないのであればコストを与えるとしているわけである。

おわりに

前回の独立型の補助スピンの活用に続き、関連型の補助スピンの活用ケースをシフト最適化の例を通じて検証した。

アニーリングの基本編はこれで終了として、今後は「より実用的な問題の検証」、「古典的な手法、アニーリング、量子ゲートマシンのTTSによる比較」を行う。

積み残し : 巡回セールスマン問題での分割法適用、時間制限付き巡回セールスマン問題