量子コンピュータ素人がIBM Quantum Challenge Fall 2021の難しかった問題を解説してみた

$$
\def\bra#1{\mathinner{\left\langle{#1}\right|}}
\def\ket#1{\mathinner{\left|{#1}\right\rangle}}
\def\braket#1#2{\mathinner{\left\langle{#1}\middle|#2\right\rangle}}
$$

この記事について

2か月ほど前にIBM Quantum Challenge Fall 2021 に参加して、せっかくなので当時解いてて難しかったな～という問題を復習がてら記事にしようかなと思います。
具体的には、IBM Quantum Challenge Fall 2021のChallenge 4 （日本語版はこちら）で出題された「Challenge 4c: 電池による収益の最適化を断熱量子計算で解く」を解説していこうと思います。
ちなみに最初に断っておきますが、私はこの問題をコンテスト中に解けませんでした涙
とはいえ、しょぼい解説をしても意味がないので、コンテスト後に公開されていた LukasBotschさんのソースコードを適宜お借りしながら、解説をしていこうと思います。
まだまだ、素人で分からないところだらけなので、もしご助言いただけたらコメントでしていただけると幸いです。

量子コンピュータ関係の他の記事は、下記で紹介しています。

概要

IBM Quantum Challenge Fall 2021 主に4つの問題が出題されました。Challenge 1～4まであるのですが、ポートフォリオ最適化であったり、VQEであったり、比較的実用面も意識した問題構成になっていました。
今回はその中でもChallenge 4 で出題された「Challenge 4c: 電池による収益の最適化を断熱量子計算で解く」を解説します。
これは電池の集積の最適化を行う問題に対して価格予測に基づいた市場の収益と予想される電池の経年劣化の両方を考慮して、各時間における電池の市場を選択するような問題をQAOAで解くといった問題になります。
面白いなあと思ったのはQAOAは通常、QUBOを使って定式化を行い、断熱量子計算を回路上でシミュレートすることで解を獲得する（参考）のですが、QUBO化を行うときに気にしないといけないのが、解きたい最適化問題の制約です。後ほど説明をしますが、QUBOはいわゆる「制約なし」の最適化問題なので、どうしても解きたい問題に制約を入れ込みたいときは制約項を加える必要があります。等式制約は単純に二乗和が最小になるような定式化を行えば良いですが、不等式制約は変数の増加が課題となります。（参考）
それをどう解決するかがポイントとなっているのがこのChallenge 4の問題でした。
概要はここまでにして実際の解説に入っていきます。

問題：電池による収益の最適化問題

近年、再生可能エネルギーの普及で．．．と説明を書こうと思いましたが、話が長くなりそうなので省略笑
細かい背景などは、こちら の「電池による収益の最適化問題」を見てみてください。
ここでは問題の定義のみを述べていきます。（こちらも上記を引用させていただきます。）

2つの市場$M_{1}$ , $M_{2}$ を考える。$n$ 個の時間枠が存在し、電池は各時間枠において、どちらか一方の市場で動作する。毎日が独立しているとみなされ、1日のうちの最適化は別の問題とみなす。充電は考慮しない。
$n$個の時間枠で2つの市場が利用可能と予測するため、各時間枠$t$（日）と各市場において、以下が既知であると仮定する：
・ 日々の収益 $\lambda_{1}^{t}$ , $\lambda_{2}^{t}$
・ 電池の日々の劣化、または健康コスト（サイクル数）$c_{1}^{t}$, $c_{2}^{t}$
上記仮定のもとで、$C_{max}$サイクル以下のコストで寿命時間と収益を最適化する。
決定変数$z_{t}, \forall t \in [1, n]$を導入し、すべての取りうるベクトル、つまりスケジュール$z = [z_{1}, ..., z_{n}]$に関して、$M_{1}$が選択された場合は$z_{t} = 0$、$M_{2}$が選択された場合は$z_{t} = 1$としたとき、

\begin{equation}
\underset{z \in \left\{0,1\right\}^{n}}{max} \displaystyle\sum_{t=1}^{n}(1-z_{t})\lambda_{1}^{t}+z_{t}\lambda_{2}^{t}
\end{equation}

\begin{equation}
    s.t. \sum_{t=1}^{n}[(1-z_{t})c_{1}^{t}+z_{t}c_{2}^{t}]\leq C_{max}
\end{equation}

とあらわすことができる。

という問題です。
目的関数だけを考えた場合、QAOAで回路を作るときは以下のような回路で表現すれば良いのですが、制約項が無視されてしまします。

参考: Knapsack Problem variants of QAOA for battery revenue optimisation

量子アニーリングですが、制約項の表現の仕方は簡単な方法でいうと以下の記事で紹介させていただいた通り、二乗和が最小になるように制約項を追加してあげてペナルティをかけることで目的関数を定義すれば良いということになります。

マクドナルドの記事を量子コンピュータでやってみた

が、これは等式制約にのみ有効な手段で、上記記事でも簡単に述べていますが、不等式制約を扱うには多くの変数が必要となるため、目的関数内に直接制約を埋め込むことができなさそうです。

ということで今回はまた別の解き方を行っていくわけになります。

Knapsack Problem variants of QAOA for battery revenue optimisation

では、実際にどうやって解いていくかを示していきます。
とはいっても実際にはあまり難しくなくて（解けなかったお前が何を言ってるんだ！というのは置いておいて笑）、章題にもある通り、実はこちらのKnapsack Problem variants of QAOA for battery revenue optimisation という論文に簡単な解き方が書いてありました。
今回のQuantum Challenge ではこれを丁寧に実装できるかがポイントだったようです。
とはいえ、これを丁寧に実装するのがまた難しいところだったので、この論文の解説とソースコード上の実装を行き来しながら、解説をしていこうと思います。

アルゴリズムの流れ

全体の流れ

アルゴリズムの全体だけ見ると通常のQAOAと同じで以下のように位相演算子 $U(C,\gamma_i)$ とミキシング演算子 $U(B,\beta_i)$ の$p$回の繰り返し構造をとります。

上記の中で大きく変わるのは 位相演算子 $U(C,\gamma_i)$の中身です。ということで、中身を紐解いていきましょう。

位相演算子とミキシング演算子

全体の流れがオリジナルのQAOAと同じことは示されたので少しずつ中身に入っていきます。
位相演算子とミキシング演算子のそれぞれでは以下の通りの手順に詳細化されます。

1. 位相演算子$U(C, \gamma_i)$
   1. 収益計算：Phase Return
   2. ペナルティー計算
      1. 劣化コスト計算：Cost Calculation
      2. 制約テスト($C_{max}$を超えたデータインデックスのマーク付け)：Constraint Testing
      3. ペナルティーの追加：Penalty Dephasing
      4. 制約テストと劣化コスト計算の再初期化(data registerとflag registerのクリーンアップ)：Initialization
2. ミキシング演算子$U(B, \beta_i)$：Mixising Operator

図で示すと以下の通りに対応しています。

完結に言うと位相演算子の中でうまくペナルティ計算をする仕組みを整えることでQAOAの枠組みの中で制約を加味した計算を行っていきます。

問題の再定義

ここで改めて、問題を定義しなおします。と言っても定式化をし直すだけで意味は同じです。
先ほど、本問題は以下のように定義されることを示しました。

\begin{equation}
\underset{z \in \left\{0,1\right\}^{n}}{max} \displaystyle\sum_{t=1}^{n}(1-z_{t})\lambda_{1}^{t}+z_{t}\lambda_{2}^{t}
\end{equation}

\begin{equation}
    s.t. \sum_{t=1}^{n}[(1-z_{t})c_{1}^{t}+z_{t}c_{2}^{t}]\leq C_{max}
\end{equation}

このままでは古典コンピュータだと解釈可能ですが、量子回路上に落とし込むことができません。
なので以下の通りに定義しなおします。

\begin{equation*}
\text{maximize } f(z)=return(z)+penalty(z)
\end{equation*}

\begin{equation*}
\text{where} \quad return(z)=\sum_{t=1}^{n} return_{t}(z) \quad \text{with} \quad return_{t}(z) \equiv\left(1-z_{t}\right) \lambda_{1}^{t}+z_{t} \lambda_{2}^{t}
\end{equation*}

\begin{equation*}
\quad \quad \quad \quad \quad \quad penalty(z)=\left\{\begin{array}{ll}
0 & \text{if}\quad  cost(z)<C_{\max } \\
-\alpha\left(cost(z)-C_{\max }\right) & \text{if} \quad cost(z) \geq C_{\max }, \alpha>0 \quad \text{constant}
\end{array}\right.
\end{equation*}

簡単に概要を述べると、

• (A) 劣化コストが$C_{max}$より小さい場合は、制約が満たされているのでそのまま目的関数を最大化します。
• (B) 劣化コストが$C_{max}$より大きい場合は、制約が超えた分だけ目的関数にペナルティをかける。

ということを行います。
先ほど、位相演算子とミキシング演算子の説明部分で示した「1. 収益計算：Phase Return」に関しては(A)においても(B)においても計算しますが、
「2. ペナルティー計算 -> 3. ペナルティーの追加：Penalty Dephasing」以降に関しては、それより手前の「2. ペナルティー計算 -> 2. 制約テスト($C_{max}$を超えたデータインデックスのマーク付け)：Constraint Testing」の評価によって、(B)では計算するし、(A)では計算しないといった制御をうまく行うことで上記のようなペナルティのかけ方を実現します。

実装

それでは前提がそろってきたので解説しながら、実装を行っていこうかと思います。
実装に関しては、LukasBotschさんのソースコード を参考に行っています。

インポート

とりあえずインポート分をまず張っておきます。(いらないものも入っています。)

from typing import List, Union
import math
from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister, assemble
from qiskit.compiler import transpile
from qiskit.circuit import Gate
from qiskit.circuit.library.standard_gates import *
from qiskit.circuit.library import QFT
simulator = Aer.get_backend('aer_simulator')

収益計算：Phase Return

では収益計算パートからいきます。時間枠、つまり、$n$が2の時を例に回路図で示すと以下の橙色の部分が該当します。

数式上は以下の箇所の評価が行われている部分となります。

\begin{equation*}
\text{where} \quad return(z)=\sum_{t=1}^{n} return_{t}(z) \quad \text{with} \quad return_{t}(z) \equiv\left(1-z_{t}\right) \lambda_{1}^{t}+z_{t} \lambda_{2}^{t}
\end{equation*}

こちらは実際のQuantum Challengeの解説にもあった通り、通常のQAOAと同様に考え、以下のように$e$の肩にそのまま評価式をのせるイメージで考えれば良いかと思います。

\begin{equation*}
\begin{aligned}
e^{-i \gamma_i . return(z)}\left|z\right\rangle 
&=\prod_{t=1}^{n} e^{-i \gamma_i return_{t}(z)}\left|z\right\rangle \\
&=e^{i \theta} \bigotimes_{t=1}^{n} e^{-i \gamma_i z_{t}\left(\lambda_{2}^{t}-\lambda_{1}^{t}\right)}\left|z_{t}\right\rangle \\
\text{with}\quad \theta &=\sum_{t=1}^{n} \lambda_{1}^{t}\quad \text{constant}
\end{aligned}
\end{equation*}

これは、それぞれの量子ビットに対して、回転ゲート $U_1\left(\gamma_i \left(\lambda_{2}^{t}-\lambda_{1}^{t}\right)\right)= e^{-i \frac{\gamma_i \left(\lambda_{2}^{t}-\lambda_{1}^{t}\right)} 2}$ で実装できます。
ということで、ソースコード上は以下の通りになります。（ここで詰まる方がもしいればこちらが参考になるかもしれません。）

def phase_return(index_qubits: int, gamma: float, L1: list, L2: list, to_gate=True) -> Union[Gate, QuantumCircuit]:
    qr_index = QuantumRegister(index_qubits, "index")
    qc = QuantumCircuit(qr_index)
    
    ##############################
    ### U_1(gamma * (lambda2 - lambda1)) for each qubit ###
    for iq in range(index_qubits):
        qc.p(-gamma * (L2[i] - L1[i])/2, iq)
    ##############################
    
    return qc.to_gate(label=" phase return ") if to_gate else qc

ペナルティー計算

ここからペナルティ計算に入ります。上記でも説明した通り、「劣化コスト計算：Cost Calculation」、「制約テスト($C_{max}$を超えたデータインデックスのマーク付け)：Constraint Testing」の計算は常に行われ、それ以降の計算の有無は劣化コストの大きさによって変わります。特にそこの分かれ目に注意しつつ、読んでいただけると幸いです。

劣化コスト計算：Cost Calculation

劣化コスト計算ではこちらの部分の計算となります。

劣化コスト計算を行うためにはまず以下の各時間枠に対するコストの足し算を計算する必要が出てきます。

\sum_{t=1}^{n}[(1-z_{t})c_{1}^{t}+z_{t}c_{2}^{t}]

ということで、まずは加算器 を定義してあげます。

加算器: QFT加算器

実は今回のQuantum Challengeでは足し算を行う上で、QFT加算器を使いましょう！と書いてありました。
ということなのでQFT加算器を使う前提の説明をしていこうかと思います。
割と何となくな説明で行くので厳密に学びたい方はこちら の参考文献を熟読いただけると幸いです。

まず、ある量子状態 $\ket{a}$ ($n+1$量子ビットを用いて$\ket{a} = \ket{0}\ket{a_1}\ket{a_2}...\ket{a_n}$と表される) をQFT（量子フーリエ変換）にかけると以下のように変換されます。

QFT\ket{a} = \frac{1}{\sqrt{2^{n+1}}}\sum_{j=0}^{2^{n+1}-1}e^{i\frac{2\pi a_k}{2^{n+1}}} \ket{k} (= \ket{\psi}) 

上記に対して逆フーリエ変換をかけるともちろん、

QFT^{-1}\ket{\psi} = \ket{a}

と元に戻ります。
$\ket{a}$ に足し算したい$\ket{b}$ ($n+1$量子ビットを用いて$\ket{b} = \ket{0}\ket{b_1}\ket{b_2}...\ket{b_n}$と表される)があるとき、仮にとある回路 $U$ として、以下のような状態を作ることができるとします。

U\ket{\psi} = \frac{1}{\sqrt{2^{n+1}}}\sum_{j=0}^{2^{n+1}-1}e^{i\frac{2\pi (a_k+b_k)}{2^{n+1}}} 

そうすると、
逆フーリエ変換をかけると、

QFT^{-1}\ket{\psi} = \ket{a + b}

という状態が作れます。つまり足し算ができちゃうわけです。そしてこれを実現する回路が存在します。
詳細はさらに参考文献になってしまうので省きますが、こちらの論文に記載の「Draper’s scheme」というものになるそうです。
ではどんな回路なのかというと以下になります。

$R_{l}$は 以下の通り定義される回路となります。

R_l = \begin{bmatrix} 1&0\\ 0&e^\frac{2\pi i}{2^l} \end{bmatrix}

また、実は今回は量子状態で足し算をしたいわけではないというところがまた一つポイントです。
上記の図は制御ゲートをうまく使って表現していますが、実はこれを表現する必要がなく、例えば、劣化コスト9を足したい。というときは、$\ket{1001} = \ket{1}\ket{0}\ket{0}\ket{1}$を足すと考え、以下の回路で済むことがわかります。

もう少し整理すると、以下の通りの固定の回路となります。

これを一旦任意の定数 const に対して行う回路は以下の通りに定義できます。

################## とあるUの定義 ####################
def subroutine_add_const(data_qubits: int, const: int, to_gate=True) -> Union[Gate, QuantumCircuit]:
        qc = QuantumCircuit(data_qubits)
        ##############################

        N = data_qubits-1

        bin_const = [int(x) for x in bin(const)[2:]]  # 向きに注意
        bin_const = [0]*(N-len(bin_const)) + bin_const

        for n in range(1, N+1):
            if bin_const[n-1]:
                qc.p(2*np.pi / 2**(n+1), N)

        for i in range(1, N+1):
            for n in range(N-i+1):
                if bin_const[n+i-1]:
                    qc.p(2*np.pi / 2**(n+1), N-i)


################ QFTとIQFTで挟む部分(上記QFT加算器の全体の定義) #####################
def const_adder(data_qubits: int, const: int, to_gate=True) -> Union[Gate, QuantumCircuit]:
    
    qr_data = QuantumRegister(data_qubits, "data")
    qc = QuantumCircuit(qr_data)
    
    ##############################
    ### QFT ###
    # Provide your code here
    for i in range(data_qubits):
        n = data_qubits - i  - 1
        qc.h(n)
        for j in range(n):
            qc.cp(pi/2**(n - j), j, n)
    ##############################

    ##############################
    ### Phase Rotation ###
    qc.add(subroutine_add_const(data_qubits, const))
    ##############################
    
    ##############################
    ### IQFT ###
    # Provide your code here
    for i in range(data_qubits):
        n = i
        for j in range(n):
            k = n - j - 1
            qc.cp(pi/2**(n - k), k, n)
        qc.h(n)
    ##############################
    return qc.to_gate(label=" [ +" + str(const) + "] ") if to_gate else qc

ちなみに上記ソースコード内にも # 向きに注意 と記載させていただいた通り、ここの演算は向きを間違えやすいです。
コンテスト本番で私はここの向きを間違えて結局最後まで直せず。。。解けなかったので皆さんは注意してみてみてください。
（ということなのでここはより正確なLukasBotschさんのソースコードを利用させていただきました。）

制御ゲートとQFT加算器を使った劣化コストの計算

QFT加算器（add(x)でxを足す回路とする）を作ることができたのでそれを使って劣化コストを足していきます。
まず、思い出してほしいのが、劣化コストは以下のように計算します。

\begin{equation}
\sum_{t=1}^{n}[(1-z_{t})c_{1}^{t}+z_{t}c_{2}^{t}]
\end{equation}

理解を簡単にするため、$z_1$ のみで考えると、

• $z_1=1$の時、$c_2^1$を加算
• $z_1=0$の時、$c_1^1$を加算

という演算を行えば良さそうです。そして量子回路では0 or 1の時、回路を適用するかどうかをうまく切り替える制御ゲートを使えば、$z_1$ のみで考えた場合、以下の通りに回路を定義すれば良いことがわかります。

注意したいのは、$z_1=0$の時、$c_1^1$を加算をしたいので $c_1^1$ 加算の前に0から1に反転させる用途でXゲートを追加し、元に戻すという用途で加算後はXゲートを再度かけなおす必要があります。
上記を踏まえ、全体の劣化コストを計算すると以下の通りになります。

これを実際に前述で定義したconst_adderを使ってコードに落とし込むと以下のようになります。

def cost_calculation(index_qubits: int, data_qubits: int, list1: list, list2: list, to_gate = True) -> Union[Gate, QuantumCircuit]:
    
    qr_index = QuantumRegister(index_qubits, "index")
    qr_data = QuantumRegister(data_qubits, "data")
    qc = QuantumCircuit(qr_index, qr_data)
    # https://qiskit.org/documentation/locale/ja_JP/tutorials/circuits_advanced/01_advanced_circuits.html

    for i, (val1, val2) in enumerate(zip(list1, list2)):
        
        ##############################
        ### Add val2 using const_adder controlled by i-th index register (set to 1) ###
        # Put your implementation here
        custom = const_adder(data_qubits, val1).to_gate().control(1)
        qc.append(custom, [qr_index[i]] + qr_data[:])
        ##############################
        qc.x(qr_index[i])
        
        ##############################
        ### Add val1 using const_adder controlled by i-th index register (set to 0)  ###
        # Put your implementation here
        custom = const_adder(data_qubits, val2).to_gate().control(1)
        qc.append(custom, [qr_index[i]] + qr_data[:])
        ##############################
        qc.x(qr_index[i])
    
    return qc.to_gate(label=" Cost Calculation ") if to_gate else qc

少し寄り道：回路の短縮

上記のcost_calculationだとcontst_addrが呼ばれるたびにQFTが計算され、回路が大きくなることがわかると思います。そのため、もう少し回路を工夫する必要があります。
とはいってもそんなに難しいことはなく、例えば、ある量子状態$\ket{a}$に$\ket{b}$だけでなく、$\ket{c}$も足したい！という場合を考えましょう。先ほど示した、$\ket{b}$の足し算を可能にするとある回路を$U_b$、$\ket{c}$の足し算を可能にするとある回路を$U_c$とすると、$\ket{a}$に$\ket{b}$と$\ket{c}$を足して、$\ket{a + b + c}$としたいときに以下のようにすれば良いことがわかります。

\ket{a + b + c} = QFT^{-1} * U_c * QFT * QFT^{-1} * U_b * QFT\ket{a} = QFT^{-1} * U_c * QFT\ket{a + b}

ただ、実際計算の順番を考えると以下のような考え方もあります。

\ket{a + b + c} = QFT^{-1} * U_c * U_b * QFT\ket{a}

間の$QFT$と$QFT^{-1}$が打ち消しあって、最初と最後に$QFT$と$QFT^{-1}$をかければ良く、間で連続して例のとある回路を適用すれば良いことがわかります。こうするとこでQFTは最初と最後にしか適用しなくて良くなるのでこのようにコードを書き換えてあげるのが回路の短縮になります。変更があるconst_adderとcost_calculationだけ示すと以下のようになります。

def const_adder(data_qubits: int, const: int, to_gate=True) -> Union[Gate, QuantumCircuit]:
    
    qr_data = QuantumRegister(data_qubits, "data")
    qc = QuantumCircuit(qr_data)
    
    ##############################
    ### Phase Rotation ###
    qc.add(subroutine_add_const(data_qubits, const))
    ##############################

    return qc.to_gate(label=" [ +" + str(const) + "] ") if to_gate else qc

def cost_calculation(index_qubits: int, data_qubits: int, list1: list, list2: list, to_gate = True) -> Union[Gate, QuantumCircuit]:
    
    qr_index = QuantumRegister(index_qubits, "index")
    qr_data = QuantumRegister(data_qubits, "data")
    qc = QuantumCircuit(qr_index, qr_data)

    ##############################
    ### QFT ###
    # Provide your code here
    for i in range(data_qubits):
        n = data_qubits - i  - 1
        qc.h(n)
        for j in range(n):
            qc.cp(pi/2**(n - j), j, n)
    ##############################

    for i, (val1, val2) in enumerate(zip(list1, list2)):
        
        ##############################
        ### Add val2 using const_adder controlled by i-th index register (set to 1) ###
        # Put your implementation here
        custom = const_adder(data_qubits, val1).to_gate().control(1)
        qc.append(custom, [qr_index[i]] + qr_data[:])
        ##############################
        qc.x(qr_index[i])
        
        ##############################
        ### Add val1 using const_adder controlled by i-th index register (set to 0)  ###
        # Put your implementation here
        custom = const_adder(data_qubits, val2).to_gate().control(1)
        qc.append(custom, [qr_index[i]] + qr_data[:])
        ##############################
        qc.x(qr_index[i])
    
    ##############################
    ### IQFT ###
    # Provide your code here
    for i in range(data_qubits):
        n = i
        for j in range(n):
            k = n - j - 1
            qc.cp(pi/2**(n - k), k, n)
        qc.h(n)
    ##############################
    
    return qc.to_gate(label=" Cost Calculation ") if to_gate else qc

上記の通り実装することで回路の節約ができました。

制約テスト(制約を超えたデータインデックスのマーク付け)：Constraint Testing

制約テストは以下の部分になります。実際に劣化コストが足されたものを評価するので、data registerと評価した結果を保存しておくflag registerが関わってきます。

ここに関しては手順を追って説明していきます。

1. まず初めに、制約を無視してとり得る劣化コストの中で最大の値を求めておきます。それを$C_{ub}$とします。（ubはupper boundです。）
2. $2^{d - 2} \leq C_{ub} < 2^{d - 1} $ となる$d$を求めます。（このdの値がdata registerの量子ビット数となります。）
3. ここで、一旦簡単に考えるため、$C_{ub} > C_{max} $で考えます。（そうでないときはいかなる$z$の組み合わせでも常に成り立つためそもそも考えなくても良いのですが、一応。）
4. $C_{diff} = 2^{d - 1} - C_{max} - 1 (> 0 \text{条件より常に成立する})$ としたとき、以下の通り制約式を変換することができます。

\begin{equation}
    \sum_{t=1}^{n}[(1-z_{t})c_{1}^{t}+z_{t}c_{2}^{t}]\leq C_{max} \iff \sum_{t=1}^{n}[(1-z_{t})c_{1}^{t}+z_{t}c_{2}^{t}] + C_{diff} \leq C_{max} + C_{diff} = 2^{d - 1} - 1 \iff \sum_{t=1}^{n}[(1-z_{t})c_{1}^{t}+z_{t}c_{2}^{t}] + C_{diff} < 2^{d - 1} 
\end{equation}

最後の変換は整数であるため、成立します。
上記制約式を整理すると、
**「① 劣化コストに② $C_{diff}$を足したものが③ $2^{d-1}$より小さい。」という制約に置き換わります。
これは$0$から $d-1$ 番目までの$d$個のdata registerを用意した時、$d-1$番目、つまり、「data registerの最上位ビットが1になると制約が違反」**と読み替えることができます。
①劣化コストは前述した通り、cost_calculationで算出でき、② $C_{diff}$は劣化コストにもう一度だけconst_adderで$C_{diff}$を足せば良いことがわかります。これらを足した結果、最上位ビットが1になれば、flag registerを1にし、制約違反とみなし、最上位ビットが0のままだとflag registerは0のままにしておくとうまく制御できます。
これらの制御はCNOTゲートを用いて実現可能なため、最終的には以下のようなソースコードとなります。

    def constraint_testing(data_qubits: int, C_max: int, to_gate = True) -> Union[Gate, QuantumCircuit]:

        qr_data = QuantumRegister(data_qubits, "data")
        qr_f = QuantumRegister(1, "flag")
        qc = QuantumCircuit(qr_data, qr_f)

        d = data_qubits
        w = 2**(d-1) - C_max - 1 # C_diffを求める
        if w > 0: # C_diffが0より大きい時、制約違反になるのでその時はflagをオンにする回路を取り入れる
            # QFT加算器でC_diffを足しこむ
            qc.append(qft(d), qr_data)
            qc.append(const_adder(d, w), qr_data)
            qc.append(qft(d, inverse=True), qr_data)
            
            # 最上位ビットが1だと制約違反なのでflag をオンにする
            qc.cx(qr_data[c], qr_f)

        return qc.to_gate(label=" Constraint Testing ") if to_gate else qc

以上が制約テストの回路となります。ここまでで制約違反を検知することができました。後はflag に基づいてペナルティ項を付け加えるだけです。

ペナルティーの追加：Penalty Dephasing

次はペナルティの追加です。以下の部分の回路です。

ここでは直前までで計算して合計劣化コストが$C_{max}$ を超えると評価された場合にペナルティを追加する回路を実装します。
説明に入る前に以下の回路を思い出してください。

\begin{equation*}
\quad \quad \quad \quad \quad \quad penalty(z)=\left\{\begin{array}{ll}
0 & \text{if}\quad  cost(z)<C_{\max } \\
-\alpha\left(cost(z)-C_{\max }\right) & \text{if} \quad cost(z) \geq C_{\max }, \alpha>0 \quad \text{constant}
\end{array}\right.
\end{equation*}

このペナルティは量子演算子$e^{i \gamma \alpha\left(cost(z)-C_{\max }\right)}$と定義できます。
ここで、$A_1[j]$を量子レジスタ$A_1$の$j$番目の量子ビット、「制約テスト(制約を超えたデータインデックスのマーク付け)：Constraint Testing」で説明した$d$を用いると、劣化コストのペナルティは以下の通り表されます。

\begin{equation*}
\alpha\left(cost(z)-C_{m a x}\right) = \alpha\left[(cost(z) + C_{diff})-(C_{max} + C_{diff}) \right] = \sum_{j=0}^{d-1} 2^{j} \alpha A_{1}[j]-(2^{d - 1} -1) \alpha
\end{equation*}

これを量子回路として表現すると、以下のようになります。

flag が1のときのみ、適用したいので、制御ゲートを用いて以下のように表現します。

ソースコード上で示すと以下のようになります。

def penalty_dephasing(data_qubits: int, alpha: float, gamma: float, to_gate = True) -> Union[Gate, QuantumCircuit]:

        qr_data = QuantumRegister(data_qubits, "data")
        qr_f = QuantumRegister(1, "flag")
        qc = QuantumCircuit(qr_data, qr_f)

        ##############################
        ### Phase Rotation ###
        # Put your implementation here
        for i in range(data_qubits):
            qc.cp((2**i) * alpha * gamma, qr_f, qr_data[i])
        
        C_max = 2**(data_qubits-1)-1
        qc.p(-C_max * alpha * gamma, qr_f)
        ##############################

制約テストと劣化コスト計算の再初期化(data registerとflag registerのクリーンアップ)：Initialization

いよいよ最後です。とはいえ、簡単な部分です。以下の箇所です。

今回の回路はQAOAなので収益計算と（違反が起きた場合、）劣化コストペナルティの和をハミルトニアンとして断熱量子計算を行っているものとみなすことができます。逆に言うとそれ以外は断熱量子計算の枠組みから取り払う、つまり、初期状態に戻す必要があります。
収益計算は phase_return、劣化コストペナルティはpenalty_dephasingで定義されているので、それ以外のcost_calculationとconstraint_testingの逆演算をかけてあげます。

   def reinitialization(index_qubits: int, data_qubits: int, C1: list, C2: list, C_max: int, to_gate = True) -> Union[Gate, QuantumCircuit]:

        qr_index = QuantumRegister(index_qubits, "index")
        qr_data = QuantumRegister(data_qubits, "data")
        qr_f = QuantumRegister(1, "flag")
        qc = QuantumCircuit(qr_index, qr_data, qr_f)
        ##############################
        
        qc.append(constraint_testing(data_qubits, C_max, False).inverse().to_gate(label="  Constraint Testing +  "), qr_data[:] + qr_f[:])
        qc.append(cost_calculation(index_qubits, data_qubits, C1, C2, False).inverse().to_gate(label="  Cost Calculation +  "), qr_index[:] + qr_data[:])

        ##############################

        return qc.to_gate(label=" Reinitialization ") if to_gate else qc

ミキシング演算子：Mixising Operator

本当の最後です。

とはいえ、QAOAで一般的に使われるミキシング演算子なのでソースコードのみで省略します。

    def mixing_operator(index_qubits: int, beta: float, to_gate = True) -> Union[Gate, QuantumCircuit]:

        qr_index = QuantumRegister(index_qubits, "index")
        qc = QuantumCircuit(qr_index)

        ##############################
        ### Mixing Operator ###
        # Put your implementation here
        qc.rx(2*beta, qr_index)
        ##############################

        return qc.to_gate(label=" Mixing Operator ") if to_gate else qc

以上で実装は終わりです。
それぞれの関数の呼出し部分などは問題としてすでに用意されていたので特に解説はしません。（ただ繰り返しで呼ぶだけです。）
興味のある方は、是非 Challenge-4 （日本語版）で動かしてみてください。

最後に

まだまだ素人でこういった応用的な問題を解くには知識が足らなすぎますが、浅い知識の中、わかる範囲で解説を書いてみました。
もし有識者の方がいれば是非、コメントいただければ泣いて喜びます。
最後まで読んでいただきありがとうございます。