トヨタ自動車 実課題プログラミングコンテスト 2023 Spring 14位解法

初めに

2023年3月5日から3月19日に開催されていたトヨタ自動車 実課題プログラミングコンテスト 2023 Springに参加し、14位になりました。

本記事では、このコンテストに対する私の解法を説明すると共に、より一般的な詰め込み問題に対する知見を整理し、また、本課題に対して更なる考察を行います。

ご参考になれば幸いです。

• トヨタ自動車 実課題プログラミングコンテスト 2023 Spring 14位解法
  • 初めに
  • 問題概要
  • 問題の性質
    • 荷物の個数
    • コンテナの体積に対する荷物の総体積の割合
  • 第1部 私の解法
    • Step1 荷物の配置位置の最適化
      • BL安定点
      • 木探索
        • ビームスタックサーチ
        • chokudaiサーチ
        • ビームサーチ
    • Step2 荷物の詰め込み順序の最適化
      • 転倒数最小トポロジカルソートの解法
      • 転倒数最小トポロジカルソートのコード
      • 実験概要
      • 実験結果
    • 解法まとめ
    • 結果
    • 考察
  • 第2部 別アプローチに関する検討
    • NFP
    • 3次元のレクトリニア図形に対するBL点
    • NFPのコード
    • 別アプローチについてのまとめ
  • 第3部 本課題に対する更なる考察
    • 定式化の問題点
      • 問題点1 荷物の詰め込みやすさ
      • 問題点2 荷物の積み重ね制限の有効性
    • 定式化の問題点についてのまとめ
  • おまけ
  • 最後に

問題概要

以下に問題概要を記します。

• コンテナに荷物を詰め込む。その向き、配置、順序を最適化せよ。
• コンテナの形状は下図に示す通り。$W \times H \times D$ の直方体の四隅に、$B \times B \times D$ の直方体ブロックが置かれた形状をしている。
• 同一種類の荷物が複数個存在する場合がある。
• 各荷物は各軸について90度単位で回転可能であるが、一部の種類の荷物は底面を固定した回転しか許されない。
• 荷物の種類によって、その上に他の荷物を重ねて置けるかどうかが制限されている。
• 荷物は、真上から垂直に下ろすことで1つずつ積み込んでいく。 このため、荷物を配置する予定位置の上部の空間には他の既に配置した荷物が存在してはならない。
• 積み込まれた荷物は、底面の面積の6割以上が、コンテナの底、もしくは他の荷物に接している必要がある。

(画像は問題文より)

また、スコア計算は、以下のペナルティによって行われます。ペナルティは低ければ低い程良いです。

(ペナルティ)=1000+\mathrm{MaxHeight}+(順番前後ペナルティ)+(収容失敗ペナルティ)

ただし、

• MaxHeightは積み込んだ荷物の上面の高さに関する最大値
• 順番前後ペナルティは荷物のindexに関する転倒数の1000倍の値
• 収容失敗ペナルティは上面の高さが $D$ を越えた荷物の体積の1000倍にほぼ比例する値

となっています。

このことから、本問はほぼ厳密に順序付けされた多目的最適化問題となっており、MaxHeight<順番前後ペナルティ<収容失敗ペナルティという順序になっています。

より詳細には、下記コンテストサイトからご覧ください。

https://atcoder.jp/contests/toyota-hc-2023spring

問題の性質

まず、入力データの性質を概観します。

0から9999までのseedの入力に関する統計情報としては、以下の通りになっていました。

荷物の個数

99%以上の確率で、荷物の個数は64個以下のようです。これにより、一部bit演算による高速化などが可能です。

コンテナの体積に対する荷物の総体積の割合

荷物の総体積の割合は、80%以下のものが大半となっています。

なお、この割合に応じて解法を少し変えることも試しましたが、私の場合、optunaによってそれは却下されました。

第1部 私の解法

それでは、私の解法を説明します。以下が最終提出です。

https://atcoder.jp/contests/toyota-hc-2023spring/submissions/39889768

私の解法は以下の二要素からなります。

1. 荷物の配置位置の最適化
2. 荷物の詰め込み順序の最適化

この二要素それぞれについて見ていきます。

Step1 荷物の配置位置の最適化

まず、荷物の配置位置の最適化について説明します。

大まかな方針としては、ビームスタックサーチ、chokudaiサーチ、ビームサーチを足して3で割った感じ、というのが表現として最も近いと思います。

BL安定点

まず、荷物を配置する箇所の絞り込み方法について説明します。

これは他の多くの参加者の方と同じで、BL安定点をベースとして候補点を列挙しました。

BL安定点の定義は以下の通りです。

アイテムを重なりなく置ける位置の中で，左方向にも下方向にも並進させることができない位置をBL安定点という．

(出典: 株式会社NTTデータ数理システム.“BL安定点”.NTT DATA.https://www.msi.co.jp/solution/nuopt/glossary/term_adf467a563c21d8012f518718e0347411a795785.html,(最終閲覧日:2023年3月24日))

より詳細には、以下の参考文献などを参照して下さい。

https://orsj.org/wp-content/or-archives50/pdf/bul/Vol.50_05_335.pdf

https://qiita.com/hari64/items/e66fd2a281ab76a58af3#%E9%95%B7%E6%96%B9%E5%BD%A2%E8%A9%B0%E3%82%81%E8%BE%BC%E3%81%BF%E5%95%8F%E9%A1%8C

梅谷 俊治.しっかり学ぶ数理最適化 モデルからアルゴリズムまで.講談社,2020.

その上で、本問には「積み込まれた荷物は、底面の面積の6割以上が、コンテナの底、もしくは他の荷物に接している必要がある」という制約があるため、単に並進不可能な点を列挙するだけでは不十分になり得ます。他の荷物の辺と合わせるようにして荷物を配置しないと、その6割を確保できない場合もあるからです。

このようなことも考慮しながら候補点を列挙しました。

関連して、四隅にある $B \times B$ のブロックも、本問特有の事情を生んでいて面白いと感じました。

通常であれば、BL点(Bottom-Left Point)、つまり、y座標→x座標の順序で最小の点を見れば事足りますが、今回はそれだけでは不都合が生じる場合があり、LB点(Left-Bottom Point)、つまり、x座標→y座標の順序で最小の点を見なければならない時もあります。

通常のBL法では面積や体積の降順で荷物を配置していくのが定石と言われていますが、今回の場合、最初から四隅に小さな $B \times B$ のブロックが存在している為、既にこの定石が破られています。なので、このような変則的なことが発生している、とも言えます。

木探索

次に、木探索の方法について述べます。

先程述べた様に、大まかな説明としては、「ビームスタックサーチ、chokudaiサーチ、ビームサーチを足して3で割った感じ」です。

それぞれの手法について、参考文献を挙げておきます。

https://qiita.com/tsukasa__diary/items/f050291b6adad80b7998

https://chokudai.hatenablog.com/entry/2017/04/12/055515

https://qiita.com/thun-c/items/058743a25c37c87b8aa4

その上で、私の解法のどのような点が各手法に由来しているかを、簡単にではありますが説明します。

ビームスタックサーチ

まず、今回の問題は、解や状態の枝刈りが極めて広範に出来る、ということが特徴として挙げられると思います。私の解法は焼き鈍し系ではなく木探索系の解法ですので、その特徴を活かすことが出来ます。

本問のペナルティは、荷物を詰め込むに従って広義単調増加しますし、また、各ペナルティの順序付けがかなり厳格なので、枝刈りが行いやすいです。

例えば、一度でも全ての荷物を収容できる解を見つけたら、それ以降は1cmでも荷物が上部からはみ出した解の探索を直ちに終了できますし、一度でも順序入れ替えなしで荷物を収容できる解を見つけたら、それ以降は1つでも順番が前後した解の探索も直ちに打ち切れます。

Shun_PIさんがコンテスト開始当初に言及されていた、枝刈り乱択という解法も大いに参考にさせて頂きました。

そして、木探索系の手法で枝刈りが重要となると、真っ先に思い出されたのがビームスタックサーチだったので、先述のtsukasaさんの記事をもう一度読み込みなおすことが出発点となりました。

ビームスタックサーチは分枝限定法にビームを導入したものである

とは、その記事で主張されていらっしゃることですが、正にそういった解釈でも理解出来ると思われます。

このような考え方に基づき、ありとあらゆる場面で徹底的に枝刈りを行うことで探索の効率化と高速化を図りました。

chokudaiサーチ

その上で、最初はほぼchokudaiサーチ的な手法を用いていました。

細いビームを何回も打つという点や、既に見た状態を一部保持するという点は、最終的な解法においても取り入れられている考え方です。

また、荷物を置く順序も、基本は体積の降順ですが、ビームごとに少しランダムな置換を加えています。

しかし、後述する工夫等の関係で、次第にビームサーチへと寄っていきました。

ビームサーチ

以上の考えを基に、ビームサーチをします。

今回行った高速化の内、特に重要な工夫として、「(遷移ごとに)状態をコピーするのではなく毎回知りたい状態をシミュレーションする」という工夫が挙げられます。

これはrhooさんの記事において言及されている手法です。実装に相違点こそありますが、根本的なアイデアは全く同じです。

https://qiita.com/rhoo/items/f2be256cde5ad2e62dde

具体的には、特定の荷物をおく、という操作を各ノードに見立てて、それを毎回計算することで、現在の荷物の配置状況を復元します。

それにより、解の情報をコピーする回数が大幅に減り、高速化に繋がります。

また、このような状態の持ち方をすることで、必要な安定点、不要な安定点を各ノード毎に効率的に持たせることも可能になり、安定点の計算を一部省略することも出来ました。

(プロファイリングの結果、この安定点の列挙が特に時間のかかる操作だと判明していたので、この部分の高速化の寄与はそれなりに大きいと思われます)

他にも、ハッシュによって重複除去をしたり、多様性確保の為に評価関数を複数用意したりした他、(これもrhooさんの記事で言及されていることですが)状態をシミュレーションする際に分かる、何段階か前の状態が一体どのノードに対応するのか、という情報をもとに、似たような形をした状態を減らす工夫などを行いました。

また、ビーム幅に関して、幅1のサーチをたくさん回す手法(chokudaiサーチ的手法)の方が良いのかな? とも思っていたのですが、optunaでパラメータ調整した結果、ビーム幅は43が最適という結果が出ました。

無論、一概に言える事ではないと思いますが、今回は先程示したような高速化手法の恩恵を、ある程度幅のあるビームサーチの方が受けやすいという事もあり、このような結果になっているものと思われます。

以上、荷物の配置位置の最適化について説明しました。

ところで、実はこの解法において一番大事な点を、あえてまだ説明していません。

それはビームサーチの評価関数をどのように設定するか、という点です。

これは私がコンテストにおいて最も長い時間をかけて取り組んだ点であり、最もスコアに大きな影響を与えた点でもあり、そして最も得られた成果が少ない点でした。

実は、最終的な私の解は、8割程度の確率で何の指標も使わず、ただ乱数に従った優先度でソートするだけ、という、それで本当にビームサーチと呼べるのかといったような評価関数が採用されることとなりました。

自分でもここが一番悔いの残る点ではあります。

この辺りの葛藤や試行錯誤は、後述させて頂きます。

Step2 荷物の詰め込み順序の最適化

続いて、荷物の配置を決めた後に、荷物の詰め込み順序の最適化を行うことを考えます。

(なお、予め断わっておくと、この部分のスコアに対する寄与は恐らく極めて軽微であり、コンテストでの順位を上げるという観点から言えば、かなり無駄な部分です。それは相対スコア制だということや、多目的最適化問題の重み付けの仕方などが原因です。しかし、現実問題としては、それなりに意味のある操作だとは思います)

この時点では、Step1によって荷物の配置が既に定められています。なので、この配置に適合する順序を、順番前後ペナルティが出来るだけ小さくなるように決定したい、というのが今の状況です。

ここで、「荷物を配置する予定位置の上部の空間には他の既に配置した荷物が存在してはならない」という制約は、「荷物Aは荷物Bの先に詰め込まれなければならない」という有向辺の形をした制約になるので、結局、以下の問題が解ければ良いです。

$N$ 頂点 $M$ 辺の有向グラフ $G$ が与えられる。各頂点 $i$ には、$\mathrm{kind}_i$ という種類が割り当てられている。 この時、$G$ のトポロジカル順序であるものの内、$\mathrm{kind}$ に関して転倒数最小の長さ $N$ の順列 $\mathrm{ord}$ を求めよ。

例を画像に示しました。左側にある図は荷物を詰め込んだ様子を横から見たものだとお考え下さい。この結果、今回は1,2,3,5,4という順列($\mathrm{ord}$)であれば、条件を満たし、かつ、転倒数は1と最小です。

ここで、画像にもある通り、この問題は転倒数最小の順列を求める代わりに、辞書順最小の順列を求める事でも、それなりに良い解を得ることが出来ます。実際、上の例では両者が一致しています。

基本的に $\mathrm{kind}$ の小さいものが前に来れば、転倒数は小さくなると考えられるので、辞書順最小の順列を求める事には妥当性があります。

また、$\mathrm{kind}$ が全ての頂点で相異なる場合、辞書順最小の順列を求めることは比較的容易です。

通常のトポロジカルソートを行う際に、入次数が0になっている頂点集合から、都度 $\mathrm{kind}$ が最小の頂点を選ぶようにすれば良いです。(このアルゴリズムが正当なことは、順列をnext_permutationで列挙して愚直解を計算し、それと照らし合わせる事でも確認しました。また、ABC223-Dが正にこの問題らしいです)

しかし、辞書順最小の順列は、必ずしも転倒数最小の順列になるとは限りません。

この例はiaNTUさんのTweetから例を借用させて頂き、掲載の許可を頂きました。ありがとうございます。

確かにそうですね、ケースが悪いです。
このケースならどうですか？
1 → 4
1 → 5
5 → 2
5 → 3

辞書順最小にすると 1 4 5 2 3 になる（転倒数4）
ですが、1 5 2 3 4にすると転倒数3が達成できる

— iaNTU (@iaNTU_) March 7, 2023

このTweetにある通り、辞書順最小の順列は、転倒数の最小値を必ず取る訳ではありません。なので、本問題は、あまり簡単な問題ではないだろうと考えられます。

余談ですが、辞書順最小の順列について、同一の $\mathrm{kind}$ が複数ある場合は、必ずしもこの方法で辞書順最小の解が得られるとは限らない事に注意して下さい。以下の例では、2と3が同一の種類の荷物です。

先述のアルゴリズムにおける、「$\mathrm{kind}$ が最小の頂点を選ぶ」ということが実行不可能な場合がある為、そのタイブレークを適切に処理しなければ、正確な辞書順最小の順列は得られません。

転倒数最小トポロジカルソートの解法

さて、以上の考察を基に、本問題に対する効率的な解答を与えます。

先程も触れた、iaNTUさんによるTweetのリプライにかなり大きなヒントを頂きましたが、bitDPをビームサーチによって簡略化することでも本問題は解くことが出来ます。

まず、通常のトポロジカルソートと殆ど同じ流れで行う点は変わりありません。

ここで、先程は入次数が0になった頂点集合から、$\mathrm{kind}$ が最小の頂点だけを選んでいましたが、今度はこれを全て試すことを考えます。(実際に全てを試す訳ではありませんが)

これまでにどの頂点を選択したか(64頂点以下の場合は、これは64bit整数にbit毎のflagを立てることによって管理出来ます)と、現在の入次数(トポロジカルソートをする上での途中経過に相当)という2つの情報があれば、bitDPによって先述の遷移を全て網羅することが可能で、厳密解を得ることが出来ます。

しかし、それだけでは実行時間が長くなりすぎる場合があるので、現在までに確定した転倒数が小さいものから上位k個のみを保持して遷移するようにすれば、ビームサーチで本問題を解くことが出来ます。

この辺りの話は、以下の記事に関連や背景等を含めて説明されている為、ビームサーチが何かをご存じでない方などはこちらをご覧ください。

https://hakomof.hatenablog.com/entry/2018/12/06/000000

(なお、別の見方をすれば、重複除去の為のhashとして、bitDPのbitに相当する整数を使用している、とも解釈できます)

また、高速に事前計算することが出来る、(擬似的な)辞書順最小の解の結果によって状態を枝刈りすることも行っています。

個人的にはヒューリスティックの非常に良い例題のように感じて面白かったです。

転倒数最小トポロジカルソートのコード

以下にコードを示します。一部は省略していますが、後述のコードも合わせてご覧ください。

using namespace std;
#include <bits/stdc++.h>

namespace minInvTopoSort {
template <class T>
struct fenwick_tree {
    fenwick_tree() : _n(0) {}
    explicit fenwick_tree(int n) : _n(n), data(n) {}
    explicit fenwick_tree(vector<T>& As) : _n(As.size()), data(As) {}
    // A[idx]+=x
    void add(int idx, T x) {
        assert(0 <= idx && idx < _n);
        idx++;
        while (idx <= _n) {
            data[idx - 1] += x;
            idx += idx & -idx;
        }
    }
    // Σ_[l,r)
    T sum(int l, int r) const {
        assert(0 <= l && l <= r && r <= _n);
        return _sum(r) - _sum(l);
    }
    // Σ_[0,r)
    T sum(int r) const {
        assert(0 <= r && r <= _n);
        return _sum(r);
    }
    // A[idx] O(logN)
    T get(int idx) const {
        assert(0 <= idx && idx < _n);
        return sum(idx, idx + 1);
    }
    // debug
    vector<T> state() const {
        vector<T> ret(_n);
        for (int i = 0; i < _n; i++) ret[i] = get(i);
        return ret;
    }

    inline T operator[](int idx) const { return get(idx); }

   private:
    int _n;
    vector<T> data;
    T _sum(int r) const {
        T s = 0;
        while (r > 0) s += data[r - 1], r -= r & -r;
        return s;
    }
};

long long inversionNumber(const vector<int>& As, int M) {
    long long ret = 0;
    fenwick_tree<int> f(M);
    for (size_t j = 0; j < As.size(); j++) {
        ret += j - f.sum(As[j] + 1);
        f.add(As[j], 1);
    }
    return ret;
}

struct TopoState {
    vector<int> indeg, zeros, ord;
    int score;
    long long hash;

    TopoState() : score(-1), hash(-1) {}
    TopoState(const vector<int>& indeg, const vector<int>& zeros)
        : indeg(indeg), zeros(zeros), score(0), hash(0) {}
    TopoState(const vector<int>& indeg, const vector<int>& zeros,
              const vector<int>& ord, int score, long long hash)
        : indeg(indeg), zeros(zeros), ord(ord), score(score), hash(hash) {}

    TopoState makeNewState(const vector<vector<int>>& adj, const vector<int>& kinds,
                           int p, int pruningScore) const {
        int newScore = score;
        for (auto& other : kinds)
            if (kinds[p] > other) newScore++;
        for (auto& prev : ord)
            if (kinds[prev] < kinds[p]) newScore--;
        if (newScore >= pruningScore) return TopoState();

        vector<int> newIndeg = indeg, newZeros = zeros, newOrd = ord;
        auto iter = find(newZeros.begin(), newZeros.end(), p);
        assert(iter != newZeros.end());
        newZeros.erase(iter);
        newOrd.emplace_back(p);
        for (int to : adj[p])
            if (--newIndeg[to] == 0) newZeros.push_back(to);

        return TopoState(newIndeg, newZeros, newOrd, newScore,
                         hash ^ (1ll << (p % 63)));
    }

    bool operator<(const TopoState& rhs) const { return score < rhs.score; }
    bool operator>(const TopoState& rhs) const { return rhs < (*this); }
};

pair<int, vector<int>> topoSort(const vector<vector<int>>& adj,
                                const vector<int>& kinds, const int M,
                                const size_t TOPO_BEAM_WIDTH) {
    size_t N = adj.size();
    vector<int> indeg(N), zeros, ord;
    for (size_t i = 0; i < N; i++)
        for (int to : adj[i]) ++indeg[to];
    for (size_t i = 0; i < N; i++)
        if (indeg[i] == 0) zeros.push_back(i);
    vector<TopoState> nowBeam;
    nowBeam.push_back(TopoState(indeg, zeros));

    // (大まかに)辞書順最小
    while (!zeros.empty()) {
        auto iter = min_element(zeros.begin(), zeros.end(),
                                [&](int a, int b) { return kinds[a] < kinds[b]; });
        int p = *iter;
        zeros.erase(iter);
        ord.emplace_back(p);
        for (int to : adj[p])
            if (--indeg[to] == 0) zeros.push_back(to);
    }
    assert(ord.size() == N);
    vector<int> Xs(N);
    for (size_t i = 0; i < N; i++) Xs[i] = kinds[ord[i]];
    int pruningScore = inversionNumber(Xs, M);
    if (pruningScore == 0 || TOPO_BEAM_WIDTH == 0) return {pruningScore, ord};

    vector<int> bestOrd = ord;

    using TopoPQ = priority_queue<TopoState, vector<TopoState>, greater<TopoState>>;
    while (!nowBeam.empty() && !nowBeam[0].zeros.empty()) {
        TopoPQ nextBeam;
        for (const auto& state : nowBeam) {
            for (int p : state.zeros) {
                TopoState newState = state.makeNewState(adj, kinds, p, pruningScore);
                if (newState.score != -1) nextBeam.push(newState);
            }
        }
        vector<TopoState> newNowBeam;
        unordered_set<long long> seen;
        while (!nextBeam.empty() && newNowBeam.size() < TOPO_BEAM_WIDTH) {
            if (seen.find(nextBeam.top().hash) == seen.end()) {
                seen.insert(nextBeam.top().hash);
                newNowBeam.push_back(nextBeam.top());
            }
            nextBeam.pop();
        }
        swap(nowBeam, newNowBeam);
    }

    int bestScore = pruningScore;
    for (auto& state : nowBeam)
        if (bestScore > state.score) {
            bestScore = state.score;
            swap(bestOrd, state.ord);
        }
    return {bestScore, bestOrd};
}
}  // namespace minInvTopoSort

実験概要

解法の性能を検証するために、以下のような実験を行いました。

$3 \leq N \leq 50$ の範囲でランダムに問題を100回生成し、それをnext_permutationで厳密に解く場合と、(擬似的な)辞書順最小の順列で解く場合、そして、様々なビーム幅で解く場合で、得られた転倒数の平均値を比較します。