補遺: 制約つきQUBOの分解法

実問題のQUBOで難しいのは、目的関数そのものよりも制約です。車は一つの経路だけを選びます。AGVは同じ通路に同時に入らないです。商品は一つの棚に一つだけ置きます。工場の工程順序は守る。資源容量は超えないです。避難経路では混雑を避けます。こうした条件をすべて二乗ペナルティで入れると、式はすぐに密になります。

制約つきQUBOの分解法は、量子アニーラを弱く見積もるための逃げ道ではありません。現場の問題を、量子アニーラが得意な二値サンプリングの部分と、古典計算が得意な連続最適化、双対更新、実行可能性検査、局所修復の部分へ分けるための設計論です。

1. 制約の二乗が作る全結合

もっとも基本的な制約は、一つだけ選ぶ制約です。対象 $i$ の候補集合を $M_i$ とし、候補 $\mu$ を選ぶ変数を $x_{\mu i}$ とすると、

$$\sum _{\mu \in M_i}{x}_{\mu i}=1$$

です。二乗ペナルティにすると

$${(\sum _{\mu \in M_i}{x}_{\mu i}-1)}^2=2\sum _{\mu <\nu }{x}_{\mu i}{x}_{\nu i}-\sum _{\mu }{x}_{\mu i}+1$$

となります。0-1変数では $x^2=x$ を使いました。ここで候補どうしの組 $\mu <\nu$ がすべて結合します。候補数が多いだけで、制約項は密なQUBOを作ります。

$K$ 個だけ選ぶ制約なら

$${(\sum _i{x}_i-K)}^2=2\sum _{i<j}{x}_i{x}_j+(1-2K)\sum _i{x}_i+K^2$$

です。この式は、同じ制約に属するすべての変数どうしを結びます。D-Wave型の量子アニーラでは物理グラフが疎なので、論理的な全結合を載せるには埋め込みが必要になります。変数一つを複数量子ビットの鎖で表すため、論理変数数よりはるかに多くの物理量子ビットを使います。

したがって、制約処理の基本方針は二つです。小さい制約は二乗ペナルティで入れます。大きい制約、広い制約、何度も出る制約は、補助場や分解で扱います。

2. 半分ずらし: 0個と1個を同じにする

交通、AGV、避難、信号制御では、「同じ場所・時刻に二つ以上入ってはいけない」という制約がよく出ます。これは「ちょうど一つ入れ」ではありません。0個でもよいし、1個でもよいです。悪いのは2個以上です。

場所・時刻の使用数を

$$n=\sum _i{y}_i$$

とします。$n=0$ と $n=1$ を同じエネルギーにし、$n=2$ から罰したいなら、

$${(n-\frac{1}{2})}^2$$

を使います。0-1変数で展開すると

$${(\sum _i{y}_i-\frac{1}{2})}^2=2\sum _{i<j}{y}_i{y}_j+\frac{1}{4}$$

です。一次項が消えることが重要です。$n=0$ と $n=1$ に差をつけず、二つ以上が同時に選ばれる組だけを罰します。

この半分ずらしは、道路混雑の見直し、AGVの衝突回避、避難経路の容量、ジョブショップの含意制約などで現れます。制約を正確に読むことは、係数をきれいにするだけでなく、目的関数に余計なバイアスを入れないためにも大切です。

3. 不等式制約とスラック変数

容量制約は

$$F(\mathit{x})\le C$$

の形です。QUBOへ素直に入れるなら、非負のスラック変数 $s$ を導入して

$$F(\mathit{x})+s=C,s\ge 0$$

とします。$s$ を二進展開

$$s=\sum _{\ell =0}^{L-1}{2}^{\ell }{y}_{\ell }$$

で表すと、ペナルティは

$$\lambda {(F(\mathit{x})+\sum _{\ell =0}^{L-1}{2}^{\ell }{y}_{\ell }-C)}^2$$

です。この方法は分かりやすいですが、スラック変数が増えます。さらに、$F(\mathit{x})$ が多くの変数の和なら、やはり広い結合が出ます。小さい容量制約には十分使えますが、大規模問題では別の方法が欲しくなります。

4. 大関法: 制約を補助場に移す

等式制約

$$F_k(\mathit{x})=C_k$$

を二乗ペナルティで入れると、

$$f(\mathit{x})=f_0(\mathit{x})+\frac{1}{2}\sum _k{\lambda }_k{\{F_k(\mathit{x})-C_k\}}^2$$

になります。大関法は、この二乗項を補助場 ${\nu }_k$ の更新へ置き換えます。固定した $\mathit{\nu }$ に対してサンプラが見る有効ハミルトニアンは

$$H_{\mathrm{e}\mathrm{f}\mathrm{f}}(\mathit{x};\mathit{\nu })=f_0(\mathit{x})-\sum _k{\nu }_k{F}_k(\mathit{x})$$

です。$F_k$ が線形和なら、$-{\nu }_k{F}_k$ は局所磁場になります。二乗ペナルティが作っていた全結合を、補助場で置き換えることができます。

固定した $\mathit{\nu }$ のもとでサンプル平均を取り、

$$⟨F_k(\mathit{x}){⟩}_{\mathit{\nu }}\simeq \frac{1}{R}\sum _{r=1}^R{F}_k({\mathit{x}}^{(r)})$$

とします。制約値 $C_k$ に近づける更新は

$${\nu }_k^{(t+1)}={\nu }_k^{(t)}+{\eta }_t\{C_k-⟨F_k(\mathit{x}){⟩}_{{\mathit{\nu }}^{(t)}}\}$$

です。サンプルが制約値より大きい方向に偏れば、その方向を抑えるように ${\nu }_k$ が動きます。サンプルが小さすぎれば、逆向きに動きます。詳しい導出は、ハバード・ストラトノビッチ変換から 補遺: ハバード・ストラトノビッチ変換と大関法 で説明しています。

5. 不等式制約に対する大関法

不等式制約

$$F_k(\mathit{x})\le C_k$$

では、制約を満たしている限りペナルティをかける必要はありません。統計力学的には、Heaviside関数を使って

$$Z=\sum _{\mathit{x}}\exp [-\beta f_0(\mathit{x})]\prod _k\mathrm{\Theta }(C_k-F_k(\mathit{x}))$$

と書けます。$\mathrm{\Theta }$ は条件を満たす配置だけを残す関数です。この関数を積分表示し、鞍点法で評価すると、非負の双対変数 ${\mu }_k\ge 0$ が現れます。固定した $\mathit{\mu }$ に対する有効ハミルトニアンは、定数を除いて

$$H_{\mathrm{i}\mathrm{n}\mathrm{e}\mathrm{q}}(\mathit{x};\mathit{\mu })=f_0(\mathit{x})+\sum _k{\mu }_k{F}_k(\mathit{x})$$

です。容量を超えやすい制約ほど ${\mu }_k$ が大きくなり、その制約に寄与する変数を選びにくくします。更新は射影つきの形になります。

$${\mu }_k^{(t+1)}={[{\mu }_k^{(t)}+{\eta }_t\{⟨F_k(\mathit{x}){⟩}_{{\mathit{\mu }}^{(t)}}-C_k\}]}_{+}$$

$[u{]}_{+}=max(u,0)$ です。制約平均が容量を超えると ${\mu }_k$ が増えます。容量に余裕があれば ${\mu }_k$ は下がり、0より小さくはなりません。これはKKT条件の相補性に対応します。

スラック変数で容量を表す方法は、制約を厳密にQUBOへ入れる方法です。一方、不等式版の大関法は、サンプラを使いながら双対変数を更新する方法です。前者は構造が単純で、後者は大規模な容量制約で結合数を抑えやすいです。

6. 列生成: 候補解を列として育てる

列生成は、大きな組合せ最適化問題を「候補解の組合せ」として解く方法です。制約つき二次0-1問題を

$$\underset{\mathit{x}}{min}\sum _{i,j}{Q}_{ij}{x}_i{x}_j\mathrm{s}.\mathrm{t}.\sum _{i,j}{A}_{kij}{x}_i{x}_j\le b_k,x_i\in \{0,1\}$$

とします。すべての0-1解を直接調べるのではなく、いくつかの候補解 ${\mathit{x}}^p$ を用意し、その凸結合

$${\mathit{x}}^a=\sum _{p\in \mathcal{P}}{\lambda }_p{\mathit{x}}^p,\sum _{p\in \mathcal{P}}{\lambda }_p=1,{\lambda }_p\ge 0$$

を考えます。制限主問題は、現在持っている候補列だけで

$$\underset{\mathit{\lambda }}{min}\sum _{p\in \mathcal{P}}\left(\sum _{i,j}{Q}_{ij}{x}_i^p{x}_j^p\right){\lambda }_p$$

を解きます。制約も

$$\sum _{p\in \mathcal{P}}\left(\sum _{i,j}{A}_{kij}{x}_i^p{x}_j^p\right){\lambda }_p\le b_k$$

と書けます。これは ${\lambda }_p$ についての線形計画です。

制限主問題の双対変数を ${\rho }_k^{\ast }$、正規化制約に対応する双対変数を ${\pi }_0^{\ast }$ とします。新しい良い候補列を探す価格付け問題は

$$\underset{\mathit{x}\in \{0,1{\}}^N}{min}\sum _{i,j}(Q_{ij}-\sum _k{\rho }_k^{\ast }{A}_{kij})x_i{x}_j-{\pi }_0^{\ast }$$

です。これはQUBOです。つまり、古典計算で制限主問題を解き、量子アニーラで新しい候補列を探します。良い列が見つかれば追加し、また制限主問題を解きます。

手順は次のように読めます。

1. 少数の実行可能候補から始めます。
2. 線形計画で、現在の候補をどう混ぜるとよいかを解きます。
3. 双対変数を使って、まだ足りない候補をQUBOとして探します。
4. 見つかった候補を列として追加します。
5. 改善できなくなるまで繰り返します。

列生成の良さは、量子アニーラに「巨大な最終問題」を渡さない点です。渡すのは、いま足りない候補を探す価格付け問題です。サンプラが多様な低エネルギー解を返すなら、それらを複数の候補列として追加できます。

7. 実行可能解を作る列生成

列生成で得られる ${\mathit{x}}^a$ は凸結合なので、一般には0-1ではありません。制限主問題の値が良くても、そのまま現場へ出せる解とは限りません。そこで、得られた候補列から実行可能な0-1解を構成する後処理が必要になります。

基本の考え方は、連続解の重み ${\lambda }_p$ を利用して有望な列を選び、0-1制約を壊さないように丸め、局所探索で修復することです。たとえば、重みの大きい列を優先して選ぶ、制約違反が小さい列の組合せを残す、QUBOサンプラで得た候補を近傍探索の初期点にする、といった操作を組み合わせます。

この後処理は地味ですが、実問題では決定的です。最適化研究では下界や緩和値が大事ですが、現場では実行可能解が必要です。量子アニーリングを列生成の価格付けに使うときも、最後に制約を満たす解へ戻す経路を持っているかが、応用の強さを決めます。

8. Benders分解: 連続変数を副問題へ送る

Benders分解は、二値変数と連続変数が混ざる問題で使いやすいです。典型形を

$$\underset{\mathit{x},\mathit{y}}{min}{\mathit{x}}^{\mathsf{T}}C\mathit{x}+{\mathit{h}}^{\mathsf{T}}\mathit{y}\mathrm{s}.\mathrm{t}.A\mathit{x}+G\mathit{y}\le \mathit{b},\mathit{x}\in \{0,1{\}}^n,\mathit{y}\ge 0$$

とします。$\mathit{x}$ を固定すると、$\mathit{y}$ については線形計画です。その双対問題は

$$\underset{\mathit{\lambda }}{max}(\mathit{b}-A\mathit{x}{)}^{\mathsf{T}}\mathit{\lambda }\mathrm{s}.\mathrm{t}.G^{\mathsf{T}}\mathit{\lambda }\le \mathit{h},\mathit{\lambda }\ge 0$$

です。双対の極点を ${\mathit{u}}_k$、極方向を ${\mathit{r}}_j$ とすると、主問題は

$$\underset{\mathit{x},t}{min}{\mathit{x}}^{\mathsf{T}}C\mathit{x}+t$$

に、切除平面

$$(\mathit{b}-A\mathit{x}{)}^{\mathsf{T}}{\mathit{u}}_k\le t$$

と実行可能性カット

$$(\mathit{b}-A\mathit{x}{)}^{\mathsf{T}}{\mathit{r}}_j\le 0$$

を加えた形になります。$t$ は「連続副問題があとで足す最小コスト」を表します。

量子アニーラで主問題を解くには、連続変数 $t$ とカットのスラックを二進展開します。たとえば

$$t(\mathit{w})=t_{min}+\delta \sum _{\ell =0}^{L_t-1}{2}^{\ell }{w}_{\ell }$$

と置きます。カット

$$(\mathit{b}-A\mathit{x}{)}^{\mathsf{T}}{\mathit{u}}_k\le t$$

は、非負スラック $s_k$ を使って

$$t(\mathit{w})+{\mathit{u}}_k^{\mathsf{T}}A\mathit{x}+s_k-{\mathit{b}}^{\mathsf{T}}{\mathit{u}}_k=0$$

とし、二乗ペナルティにします。$s_k$ も二進展開します。結果として、主問題は

$$H_{\mathrm{m}\mathrm{a}\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{r}}={\mathit{x}}^{\mathsf{T}}C\mathit{x}+t(\mathit{w})+\sum _k{P}_k{(t(\mathit{w})+{\mathit{u}}_k^{\mathsf{T}}A\mathit{x}+s_k-{\mathit{b}}^{\mathsf{T}}{\mathit{u}}_k)}^2+\cdots$$

というQUBOになります。$\cdots$ には実行可能性カットのペナルティが入ります。

Benders分解の流れは明快です。量子アニーラで二値主問題を解きます。得られた $\mathit{x}$ を固定して古典線形計画を解きます。副問題が実行可能で最適値が分かれば最適性カットを加えます。実行不可能なら実行可能性カットを加えます。主問題と副問題の差が十分小さくなれば止めます。

$$|t-(\mathit{b}-A\mathit{x}{)}^{\mathsf{T}}{\mathit{u}}_k|\le \epsilon$$

この方法は、連続変数まで無理に量子アニーラへ押し込まないところが強いです。量子アニーラは、離散的な設計選択を担当します。連続副問題は、古典最適化の安定した道具に任せます。

9. ラグランジュ分解協調: 工程を分けて会話させる

鉄鋼製造スケジューリングのように、二つの工程が結合している問題では、工程ごとにQUBOを分けると扱いやすくなります。工程 $p=1,2$、製品 $i$、スロット $s$ に対して $x_{p,i,s}$ を用意し、時刻を

$$t_{p,i}=\sum _{s}s{x}_{p,i,s}$$

とします。工程順序制約は $t_{1,i}\le t_{2,i}$ です。これをラグランジュ乗数 ${\lambda }_i$ で緩めると、

$$E_{\mathrm{L}\mathrm{D}\mathrm{C}}=\mathrm{C}\mathrm{o}\mathrm{s}\mathrm{t}+\rho \mathrm{P}\mathrm{e}\mathrm{n}\mathrm{a}\mathrm{l}\mathrm{t}\mathrm{y}+\sum _i{\lambda }_i(t_{1,i}-t_{2,i})$$

です。これは

$$E_1={\mathrm{C}\mathrm{o}\mathrm{s}\mathrm{t}}_1+\rho {\mathrm{P}\mathrm{e}\mathrm{n}\mathrm{a}\mathrm{l}\mathrm{t}\mathrm{y}}_1+\sum _i{\lambda }_i{t}_{1,i},E_2={\mathrm{C}\mathrm{o}\mathrm{s}\mathrm{t}}_2+\rho {\mathrm{P}\mathrm{e}\mathrm{n}\mathrm{a}\mathrm{l}\mathrm{t}\mathrm{y}}_2-\sum _i{\lambda }_i{t}_{2,i}$$

へ分かれます。工程1が遅れすぎるなら ${\lambda }_i$ を増やし、工程1でその製品を早める圧力、工程2で遅める圧力をかけます。更新はたとえば

$${\lambda }_i^{\mathrm{n}\mathrm{e}\mathrm{w}}={\lambda }_i+\alpha max(t_{1,i}^{\ast }-t_{2,i}^{\ast },0)$$

です。二つの工程が、それぞれのQUBOを解きながら、ラグランジュ乗数を通して会話します。工程を分けることで、変数数はおおよそ $2N^2$ から $N^2$ 規模へ下がります。

10. クラスタ分解と候補解合成

AGV、避難、配送のような空間的問題では、地理的に近い対象をクラスタに分けることができます。クラスタごとに局所QUBOを解き、複数の候補解を得ます。その後、クラスタ間の干渉を考えながら候補解を一つずつ選びます。

クラスタ $m$ の候補解 $n$ を選ぶ変数を $s_{mn}$ とします。クラスタ内コストを $A_{mn}$、クラスタ $m$ の候補 $n$ とクラスタ $p$ の候補 $q$ を同時に選ぶときの干渉を $B_{mnpq}$ とすると、

$$H(\mathit{s})=\sum _{m,n}{A}_{mn}{s}_{mn}+\sum _{m\ne p}\sum _{n,q}{B}_{mnpq}{s}_{mn}{s}_{pq}+\mathrm{\Gamma }\sum _m{(\sum _n{s}_{mn}-1)}^2$$

です。この上位QUBOは、元問題のすべての車両や避難者を直接変数にするより小さくなります。

量子アニーラの多様なサンプルは、この方法と相性が良いです。クラスタ内で一つの最良解だけを残すと、クラスタ間で組み合わせたときに詰まることがあります。少し悪いが別の経路を使う候補を残しておくと、全体では混雑が減ることがあります。サンプラは、最終解を一つ出すだけでなく、組合せ直すための候補集合を作ります。

11. どの方法を選ぶか

制約つきQUBOでは、最初から最強の分解法を使う必要はありません。まず制約の形を見ます。

一つだけ選ぶ制約が小さいなら、二乗ペナルティで十分です。二つ以上禁止なら、半分ずらしが自然です。合計や容量が大きく、全結合が問題になるなら、大関法や不等式版の双対更新を考えます。候補解の集合を育てたいなら列生成が合います。連続変数が多いならBenders分解が合います。工程や地域のように自然な分割があるなら、ラグランジュ分解やクラスタ分解が合います。

実装では、次の流れが役に立ちます。

1. 変数を決めます。
2. 絶対制約と評価関数を分けます。
3. 小さい制約はペナルティで入れます。
4. 大きい制約は双対変数、列生成、Benders、クラスタ分解へ逃がします。
5. サンプルを後処理し、必ず実行可能性を検査します。
6. 実機やサンプラの係数範囲に合わせてスケーリングします。
7. 古典ソルバ、SA、SQA、局所探索と比較します。

量子アニーリング応用は、QUBOを書いて終わりではありません。QUBOを書き、縮め、分け、サンプルし、戻し、検査するまでが一つの方法です。その全体を設計できると、量子アニーラは現場の制約の中で使える道具になります。

不等式制約に対する劣勾配法は文献 [8]を参照してください。

参考文献

1. Masayuki Ohzeki, “Breaking limitation of quantum annealer in solving optimization problems under constraints,” Scientific Reports 10, 3126 (2020). 出典

2. Sota Hirama, Masayuki Ohzeki, “Efficient Algorithm for Binary Quadratic Problem by Column Generation and Quantum Annealing,” Journal of the Physical Society of Japan 92, 113002 (2023). 出典

3. Taisei Takabayashi, Naoki Maruyama, Takuma Yoshihara, Renichiro Haba, Masayuki Ohzeki, “Efficient Construction of Feasible Solutions in Column Generation Using Quantum Annealing,” Journal of the Physical Society of Japan 94, 094005 (2025). 出典

4. Takuma Yoshihara, Masayuki Ohzeki, “Accelerating Extended Benders Decomposition with Quantum-Classical Hybrid Solver,” Journal of the Physical Society of Japan 95, 033001 (2026). 出典

5. Kouki Yonaga, Masamichi Miyama, Masayuki Ohzeki, Koji Hirano, Hirokazu Kobayashi, Tetsuaki Kurokawa, “Quantum Optimization with Lagrangian Decomposition for Multiple-process Scheduling in Steel Manufacturing,” ISIJ International 62, 1874–1880 (2022). 出典

6. Thinh Nguyen Quang, Kosuke Matsuyama, Keisuke Shimizu, Hiroki Sugano, Eiji Kurimoto, Masahiro Miki, et al., “Quantum annealing-based route optimization for commercial AGV operating systems in large-scale logistics warehouses,” Scientific Reports 15, 44047 (2025). 出典

7. Reo Shikanai, Renichiro Haba, Yusuke Okazaki, Kazumichi Matsumoto, Masayuki Ohzeki, “Large-scale evacuation route optimization leveraging sampling diversity in quantum annealing,” Scientific Reports 16, 22108 (2026). 出典

8. Taisei Takabayashi, Takeru Goto, Masayuki Ohzeki, “Subgradient Method Using Quantum Annealing for Inequality-Constrained Binary Optimization Problems,” Journal of the Physical Society of Japan 94, 054003 (2025). 出典