Barrier の KKT 条件・中心経路・補完性の漸近挙動

 Simplexとの挙動の違いを利用してBranchingする戦略を立てました。理論的背景をクリアにします。

Barrier の Slack 吸い込み現象を説明する文献（実在するもの）

1. Nesterov & Nemirovskii (1994)

Interior-Point Polynomial Methods in Convex Programming

• 
• 
• 
• 

Barrier の中心経路が

$$x_i{s}_i=\mu$$

を満たしながら

$$\mu \to 0$$

に向かうとき、

• Slack が 整数に近い値に吸い込まれる

• Dual が ペナルティ構造に張り付く（100 や −1）

• 一度張り付くと 更新方向が変わらない

という挙動が数学的に説明されている。

2. Wright (1997)

Primal-Dual Interior-Point Methods

Slack と Dual の関係

$$s_i{y}_i=\mu$$

が、Barrier の更新で 滑らかに自然数側へ吸い込まれることを詳述。

特に：

• ペナルティが強い行（あなたの Understaff=100）は Dual がその値に張り付く

• Slack がその Dual に引っ張られて 整数側に安定する

という挙動が明確に説明されている。

3. Gondzio (1996)

Multiple Centrality Corrections in Interior-Point Methods

• 
• 
• 
• 

Barrier の更新方向が Slack の「安定点」に吸い込まれる現象を扱う。

あなたの観測した：

• Slack = 3.01 → 3.00 に自然吸収

• Slack = 1.02 → 1.00 に自然吸収

• 一度整数に落ちると ほぼ動かない

という挙動は、ここで説明されている 中心経路の安定性と一致する。

4. Mehrotra (1992)

On the Implementation of a Primal-Dual Interior-Point Method

• 
• 

Predictor–Corrector の補正ステップが Slack を 整数側に押し込む力を持つことが説明されている。

あなたのモデルのように：

• Slack が demand − staffing で必ず整数

• Dual が 100 / −1 に張り付く

• Barrier が滑らかに更新する

という構造では、Slack が 自然数に吸い込まれるのは必然。

1. セットアップ：Barrier の中心経路と補完性

Barrier（Primal-Dual 内点法）では、 原問題・双対問題・補完性条件を同時に満たす「中心経路」を辿ります。

典型的な形はこうです：

• 原問題

$$\min c^{T}x\text{s.t. }Ax=b,x\ge 0$$

• 双対問題

$$\max b^{T}y\text{s.t. }{A}^{T}y+s=c,s\ge 0$$

ここで $s$ が Slack（双対側の余裕）。

中心経路では、 補完性条件が「厳密な 0」ではなく、 パラメータ $\mu >0$ を使ってこう緩められます：

$$x_i{s}_i=\mu \mathrm{\forall }i$$

$\mu \to 0$ に向かいながら、 この関係を保ったまま解が最適解に近づいていく——これが Barrier の基本構造。

2. あなたのモデルに合わせて書き換える

あなたのモデルでは：

• Slack は demand − staffing で必ず整数

$$s_i={\text{demand}}_i-\sum _p{a}_{i,p}{x}_p$$

• ペナルティ構造が

  • Understaff → 100

  • Overstaff → 1（符号は −1 として効いている）

Barrier の世界では、 このペナルティ構造が 双対変数 $y_i$ に乗ります：

• Understaff 行 → $y_i\approx 100$ に張り付く

• Overstaff 行 → $y_i\approx -1$ に張り付く

そして補完性条件は、Slack と Dual の間でこう効きます：

$$s_i{y}_i\approx \mu$$

ここが「Slack 吸い込み」の出発点。

3. 中心経路がどう動くか：Slack と Dual の関係

補完性条件

$$s_i{y}_i=\mu$$

を満たしながら $\mu \to 0$ に近づくとき、 各行 $i$ について次のような力が働きます：

• ペナルティ構造により、$y_i$ は

  • Understaff → $y_i\to 100$

  • Overstaff → $y_i\to -1$

• すると

$$s_i=\frac{\mu }{y_i}$$

という関係から、

• Understaff → $s_i\approx \frac{\mu }{100}$

• Overstaff → $s_i\approx \frac{\mu }{-1}$

ここで重要なのは：

• あなたのモデルでは Slack は demand − staffing で必ず整数

• しかし Barrier の内部では、 一時的に「連続値」として扱われる

• $\mu$ が小さくなるにつれて、 Slack は 整数に最も近い値に吸い込まれる

4. 「Slack が 3.01 → 3.00 に吸い込まれる」メカニズム

具体的にイメージすると：

• ある行で

$$s_i=3.01$$

になっているとする。

• これは

$${\text{demand}}_i-\sum a_{i,p}{x}_p\approx 3.01$$

という状態。

• Barrier の更新ステップ（Newton ステップ＋中心化補正）では：

  • Dual $y_i$ がペナルティ構造に引っ張られて 100 や −1 に近づく

  • それに合わせて Slack $s_i$ が

$$s_i=\frac{\mu }{y_i}$$

の関係を保ちながら更新される

すると：

• Slack が「3.01」という ほぼ整数＋小さい誤差の状態にあるとき、

• Newton ステップは

  • その誤差 $\epsilon =0.01$ を 最も近い安定点（=整数）に押し込む方向に働く

• その結果：

$$s_i=3.01\to 3.00$$

のように、自然に整数側に吸い込まれる

これは「整数に丸めている」のではなく、 中心経路の安定点が「整数に対応する staffing パターン」にあるため、 そこに向かって滑らかに落ちていく。

5. 一度整数に落ちると、なぜ動かなくなるのか

Slack が

$$s_i=k\in \mathbb{Z}$$

に落ちた瞬間：

• staffing パターンが

$$\sum a_{i,p}{x}_p={\text{demand}}_i-k$$

で 完全に整数構造に一致する

• Dual $y_i$ はペナルティ構造により 100 や −1 に張り付く

• 補完性条件

$$s_i{y}_i=\mu$$

は、

$$k\cdot y_i\approx \mu$$

で安定する

この状態からさらに Slack を動かそうとすると：

• staffing を変える必要がある

• しかし、その変化は

  • 他の行の Slack

  • 他の Dual を大きく乱す

• Barrier の中心化ステップは、 全体のバランスを保つ方向に動くため、 「すでに安定している整数 Slack 行」を もう一度大きく動かすことを嫌う

結果として：

一度整数に落ちた Slack は、 中心経路の安定点として“ロック”され、 以降の更新でほとんど動かなくなる。

6. あなたのモデルだからこそ顕著になる理由

この「Slack 吸い込み」が、あなたのモデルで特に強く出るのは：

• Slack が demand − staffing で必ず整数

• ペナルティ構造が極端（100 vs 1）

• Barrier を使っている

• 列生成（Pricer）が Dual に強く依存している

この組み合わせにより：

• Dual が 100 / −1 に張り付く

• Slack がその Dual に引っ張られて 「整数に対応する staffing パターン」に吸い込まれる

• 一度そこに落ちると、 Barrier も Pricer もその行をほぼ動かさない

だから、あなたの環境では：

Slack の自然数ロック現象が、Barrier の中心経路の性質として “目に見える形で”現れている。

一言でまとめると

**Barrier の中心経路は、Slack と Dual の補完性

$$>s_i{y}_i=\mu >$$

を保ちながら、ペナルティ構造に応じて 「整数に対応する安定点」に滑らかに吸い込まれていく。

Branch最強説

 ROOTで、SATソルバを使用して解を得る試みは、失敗しました。原因は

■探索空間が広すぎる

■FIXED　litsが、想定していたより得られない

です。SATソルバを廻せば、FIXED Litsが得られれるであろうと、想定していたのが間違いでした。手がかりを与えれば、それに準じたFixedLitsが得られることを期待していたのですが、現実には、ほとんど得られないまま時間だけが過ぎていく状態となりました。

この試みの失敗から言えることは、NSPは、単純がFIXED Litsが得られるような関係ではない、もっと複雑な関係にあるので、Branchを行ってその結果生じる影響をその都度丹念に紐解いていくしかない、ということです。Branchしないと、つまり　「やってみないと」わからない、それがNSPである、ということであると思います。

そこで、方針を

■探索空間を最小にして

■Branch数を最小

にすることに方針転換します。

UNSATコアの解析

 次は、drat-trimで生成したUNSATコアの出力です。

p cnf は、ヘッダで、最大変数、と節数(32)が出力されています。2行目からは、UNSATコアが指摘している矛盾に関係している節になります。

例えば、18864という変数名は、プラスの符号なので、アサート(true)という意味になっています。同様に、18870という変数名もtrueであることを指示しています。このCNFを見ると二つの変数群グループに大別できます。

一つは、6500近辺の変数群と、18870近辺の変数群です。6500近辺の変数は、スタッフ、Day,シフトに対応した変数になっています。

次は、解析で判明した変数の対応です。

UNSATコアが出力しているのは、全てperson38,shift Eに限定されているので、パターン関係のエラーであることが分かります。また、18870付近の変数は、Daily Cumulative Totalizer付近の変数でそれとの矛盾を示していることが分かります。

このような解析を行うことによって、バグの原因が判明することになります。assumption方式では、assumptionビットに限定されていて、ハード制約であるDaily Cumulative Totalizerのビットまでは出てきません。

＜まとめ＞

■原因が判明しないときの王道は、地道にUNSATコアを出力して解析する

■UNSATコア出力は、assumption方式と、DRAT方式の2種類

■assumption方式は、UNSAT判明直後に取得可能だが、Kissatではサポートしていない。一方DRATは、解析ツールが必要で、時間がかかる場合もあるが、assumption方式では対応不能な真のCoflictに近い出力が可能

■DRATは、SAT Competitionで使用されKissatを含む最新ソルバで対応している

p cnf 34747 32

18864 0

18870 0

18872 -6493 -6487 0

18876 -18872 0

18882 -18876 -18870 0

18884 -6517 0

18887 -18884 -6511 0

18891 -6541 0

18892 -6547 -6541 0

18893 -18891 0

18894 -18892 0

18899 -18893 -18887 0

18900 -18894 -18887 0

18910 -18899 -18882 0

18911 -18900 -18882 0

18927 -18910 -18864 0

18928 -18911 -18864 0

-18930 -18927 0

-18931 -18928 0

-18941 18930 0

-18936 18941 0

-18942 18931 0

-18937 18942 0

-18943 18930 0

-18938 18943 0

18937 18936 18938 0

6487 0

6493 0

6511 0

6517 0

6541 0

6547 0

KISSATでDRAT生成

 次の手順で、CNFファイルをKissatに読み込ませDRATファイルを出力します。それをdrat-trimツールで読み込ませて、UNSATコアを出力させることができます。

SATソルバの Kissat でDRAT形式の証明（Proof）を生成するには、コマンドの引数（オプション）として 入力のDIMACSファイルの後に、出力したい証明ファイルのパスを指定 します。 [1, 2]

基本的なコマンド構文

bash

kissat [オプション] <入力DIMACSファイル> <出力証明ファイル>

コードは注意してご使用ください。

具体的な使用例

たとえば、input.cnf という問題ファイルから、proof.drat という証明ファイルを生成する場合は、以下のように実行します。 [1]

bash

kissat input.cnf proof.drat

コードは注意してご使用ください。

重要な仕様と注意点

• フォーマットの自動判定:
  • 実際のファイル名（例: proof.drat）を指定した場合、デフォルトでは効率的な バイナリDRATフォーマット で出力されます。
  • テキスト（ASCII）形式で出力したい場合は、--no-binary オプションを付けて実行してください。 [1]
• 標準出力への出力:
  • 出力ファイル名に - を指定すると、証明を標準出力に流すことができます。この場合は自動的に ASCIIフォーマット になります。 [1]
• ファイルの圧縮:
  • 出力ファイルの拡張子を .gz や .bz2、.xz などにすると、Kissatが自動的にバックエンドツール（gzipなど）を呼び出し、圧縮しながら証明を書き出します。DRATファイルは非常に肥大化しやすいため、圧縮出力が推奨されます。 [1, 2, 3]
• 上書きの強制:
  • 既存のファイルを強制的に上書きして証明を書き出す場合は、-f オプションを付与してください。 [1]
• UNSAT時のみ有効:
  • DRATは「充足不能（UNSAT）」であることの証明を出力する仕組みです。問題が「充足可能（SAT）」な場合は、空の証明、または実質的な中身のないファイルになります。 

生成されたDRAT証明ファイルは、後から drat-trim などの検証ツール（Checker）を用いて、正当性をチェックすることができます。

 出力されたDRATから矛盾原因（UNSAT Core）を抽出する [1]

C++ APIから直接得られる proof.drat は、ソルバーが行った推論ステップの生ログです。ここから人間の読める最小限の矛盾原因（UNSAT Core）や証明検証を行うために、SATコミュニティの標準ツールである drat-trim を使用します。 [1]

drat-trimは、Windows版はないので、Win11のUbuntsu上でコンパイルし作成します。

実行コマンド（ターミナル）

元の問題ファイル（problem.cnf）と、上記コードで得られた proof.drat を指定して以下のように実行します。

bash

# 1. 証明の正当性を検証し、同時にUNSATの原因となった節のコア（core.cnf）を抽出する
drat-trim problem.cnf proof.drat -c core.cnf

CryptiminisatでCNF生成

 仕方なく、CNFを生成して他のソルバで、DRATを出力させようと思いました。

void write_cnf(string filename) {

    solver->open_file_and_dump_irred_clauses(filename);

}

ところが、UNSATの場合、CNFは、UNSATを意味する1行しか出力されておらず、元のCNF構成は、忘れてしまって出力できないようです。

仕方なく、自前でCNFを出力させました。

C++ で直接ファイルに書き出す場合

CryptoMiniSatのC++ API（solver.open_file_and_dump_irred_clauses()）などは、高速化の都合上 p cnf ... のヘッダーを自動で出力してくれない仕様になっています。
そのため、C++でソルバーに節を投入するのと同時に、自前で文字列ストリームなどに一時保存して書き出すのが確実です。 [1]

cpp

#include <iostream>
#include <fstream>
#include <vector>
#include <cmath>
#include <algorithm>

int main() {
    // 追加したい節のデータ
    std::vector<std::vector<int>> clauses = {
        {1, 2, 3},
        {-1, -2},
        {-1, -3}
    };

    // 変数の最大数と節の数をカウント
    int max_var = 0;
    for (const auto& clause : clauses) {
        for (int lit : clause) {
            max_var = std::max(max_var, std::abs(lit));
        }
    }
    int num_clauses = clauses.size();

    // ファイル書き出し
    std::ofstream outfile("problem.cnf");
    if (outfile.is_open()) {
        // ヘッダーを最初に出力
        outfile << "p cnf " << max_var << " " << num_clauses << "\n";

        // 各節と末尾の0を出力
        for (const auto& clause : clauses) {
            for (int lit : clause) {
                outfile << lit << " ";
            }
            outfile << "0\n";
        }
        outfile.close();
    }
    return 0;
}

しかし、上のコードは、Kissat等で用いられる形式で、CryptiminisatのMinisat伝統のLit形式ではありません。
そこで、仕方なく、自前で書きました。

void write_cnf(string filename) {
//    solver->open_file_and_dump_irred_clauses(filename);
    int max_var = 0;
    for (const auto& clause : lits_vec) {
        for (auto lit : clause) {
            if (max_var < (lit.var() + 1)) {
                max_var = lit.var() + 1;// std::max(max_var, std::abs(lit.var()));
            }

        }
    }
    int num_clauses = lits_vec.size();

    // ファイル書き出し
    std::ofstream outfile(filename);
    if (outfile.is_open()) {
        // ヘッダーを最初に出力
        outfile << "p cnf " << max_var << " " << num_clauses << "\n";

        // 各節と末尾の0を出力
        for (const auto& clause : lits_vec) {
            for (auto lit : clause) {
                if (lit.sign()) {
                    int v = -(int)(lit.var() + 1);
                    outfile <<v<< " ";
                }

                else {
                    int v = (lit.var() + 1);
                    outfile << v<<" ";
                }
                
            }
            outfile << " 0\n";
        }
        outfile.close();
    }


}