How to program a MIP solver to find balanced Gray code for mixed radices?

Question

The permutations of a mixed radix number can be ordered to achieve Grayness (in the sense of Gray code) with optimal balance and span length. Each of these constraints will be explained in turn. In my examples, I use a mixed radix number consisting of a base 2 digit, a base 3 digit, and a base 4 digit. This set is called [234], and it has 2 × 3 × 4 = 24 permutations. The permutations are listed below, in ascending order. For compactness, the digits are shown as rows, with the top row corresponding to the set’s first digit. The leftmost column is the first permutation 000, the next column is the second permutation 001, then 002, 003, 010, 011, 012, 013, and so on.

2:  0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1
3:  0 0 0 0 1 1 1 1 2 2 2 2 0 0 0 0 1 1 1 1 2 2 2 2
4:  0 1 2 3 0 1 2 3 0 1 2 3 0 1 2 3 0 1 2 3 0 1 2 3

In the above set, multiple digits may change from one permutation to the next. For example, between the fourth and fifth permutations (003 and 010), two digits change at once. To make a Gray set, we must reorder the permutations so that only one digit changes at a time. This constraint includes the wraparound from the first to the last permutation. Below is [234] reordered to be Gray:

2:  0 1 1 0 0 1 1 0 0 1 1 1 0 0 0 1 1 1 0 0 1 1 0 0 
3:  0 0 1 1 2 2 2 2 0 0 1 1 1 1 1 1 0 2 2 2 2 0 0 0 
4:  0 0 0 0 0 0 1 1 1 1 1 2 2 1 3 3 3 3 3 2 2 2 2 3

The above set is Gray, but not balanced. To be balanced, each of the set’s digits must change the same number of times, or as close as possible. In the above set, the 2’s place changes 10 times, the 3’s place changes 7 times, and the 4’s place also changes 7 times. A set’s imbalance is the absolute value of the difference between its minimum and maximum digit changes, in this case 10 – 7 = 3. Below is [234] reordered to have optimal balance; each digit changes 8 times, so the imbalance is now zero:

2:  0 1 1 0 0 1 1 0 0 1 1 1 0 0 0 1 1 1 1 1 0 0 0 0 
3:  0 0 1 1 2 2 2 2 0 0 1 1 1 1 1 1 2 0 0 2 2 2 0 0 
4:  0 0 0 0 0 0 1 1 1 1 1 2 2 1 3 3 3 3 2 2 2 3 3 2

The above set is Gray and balanced, but digits get stuck for longer than we’d like. For example, the 4’s place stays zero for the first six permutations. This constitutes a span, with a length of six. In the above set, the maximum span length is six. For optimal granularity, the maximum span should be as short as possible. Below is [234] reordered so that the maximum span length is four instead of six:

2:  0 1 1 0 0 1 1 0 0 1 1 1 1 0 0 0 0 1 1 1 1 0 0 0 
3:  0 0 1 1 1 1 0 0 0 0 1 2 2 2 2 1 1 1 0 2 2 2 2 0 
4:  0 0 0 0 1 1 1 1 2 2 2 2 0 0 2 2 3 3 3 3 1 1 3 3

The above ordering of [234] is Gray, optimally balanced, and minimizes stuck digits. It’s the best we can do for this particular set. But for larger sets, such as [345], optimal solutions are much harder to find, because my code is too slow. Can a MIP solver do better? The solution should be coded in a language that’s supported by one of the solvers available at NEOS, because these are the only high-quality solvers I have access to (for example CPLEX via AMPL, GAMS, LP, MPS, or NL). The application is atonal music theory, hence only sets with ranges of twelve or less are relevant. The complete list of sets I’m trying to optimize is here.

EDIT: Some commenters asked about my code, so I'm enclosing it below. I use Visual Studio 2012, but this code should compile fairly easily in any C++ compiler. I use x64 (64-bit code).

// Copyleft 2023 Chris Korda
// This program is free software; you can redistribute it and/or modify it
// under the terms of the GNU General Public License as published by the Free
// Software Foundation; either version 2 of the License, or any later version.

// BalaGray.cpp : Defines the entry point for the console application.
// This app computes balanced Gray code sequences, for use in music theory.

#include "stdafx.h" // precompiled header
#include "stdint.h" // standard sizes
#include "vector"   // growable array
#include "fstream"  // file I/O
#include "assert.h" // debugging

using namespace std;

#define MORE_PLACES 0   // set non-zero to use more than four places
#define DO_PRUNING 1    // set non-zero to do branch pruning and reduce runtime

class CBalaGray {
public:
// Construction
    CBalaGray();

// Attributes
    int     GetPermCount() const { return static_cast<int>(m_arrPerm.size()); }

// Operations
    void    Reset();
    void    Calc(int nPlaces, const uint8_t *arrRange);

protected:
// Constants
    enum {
#if MORE_PLACES
        MAX_PLACES = 8,
#else
        MAX_PLACES = 4,
#endif
        MAX_RANGE = 255,
        ULONGLONG_BITS = 64,
    };
    enum {  // pruning thresholds may require manual tuning; see notes in set list
        PRUNE_MAXTRANS = 18,
        PRUNE_IMBALANCE = 3,    
    };

// Types
    union PERM {    // permutation; size depends on MAX_PLACES
        uint8_t b[MAX_PLACES];  // array of places
#if MORE_PLACES
        uint64_t    dw; // double word containing all places
#else
        uint32_t    dw; // double word containing all places
#endif
    };
    struct STATE {  // crawler stack element
        uint8_t iPerm;      // permutation index
        uint8_t iGray;      // Gray neighbor index
        PERM    nTrans;     // transition counts, one per place
    };
    typedef vector<PERM> CPermArray;
    typedef vector<STATE> CStateArray;
    typedef vector<uint8_t> CPlaceArray;    // enough for atonal music theory

// Member data
    int     m_nPlaces;  // number of places
    int     m_nGrayPerms;   // number of Gray permutations reachable from a permutation
    int     m_nGrayStrideShift; // stride of Gray permutations array, as a shift in bits
    CPlaceArray m_arrRange; // array of ranges, one for each place
    CPermArray  m_arrPerm;  // array of permutations
    CPlaceArray m_arrGray;  // 2D table of permutations reachable from each permutation
    CStateArray m_arrState; // array of states; crawler stack
    ofstream    m_fOut; // output file

// Helpers
    int     Pack(const PERM& perm) const;
    void    MakePerms(int nPlaces, const uint8_t *arrRange);
    void    MakeGrayTable();
    void    DumpGrayTablePerms() const;
    void    DumpPerm(const PERM& perm) const;
    void    DumpPerms() const;
    void    DumpSet() const;
    void    WriteBalanceToLog(int nImbalance, int nMaxTrans, int nMaxSpan);
    void    WriteSequenceToLog(int iDepth);
    bool    IsGray(PERM p1, PERM p2) const;
    int     ComputeBalance(int iDepth, int& nMaxTrans, PERM& nTransCounts) const;
    int     ComputeMaxSpan(int iDepth) const;
};

CBalaGray::CBalaGray()
{
    m_fOut.open("BalaGrayIter.txt", ios_base::out); // open output file
    assert(m_fOut != NULL);
    Reset();
}

void CBalaGray::Reset()
{
    m_nPlaces = 0;
    m_arrRange.clear();
    m_arrState.clear();
}

int CBalaGray::Pack(const PERM& perm) const
{
    int nPacked = perm.b[m_nPlaces - 1];    // init total to first place
    for (int iPlace = m_nPlaces - 2; iPlace >= 0; iPlace--) {   // for each subsequent place
        nPacked *= m_arrRange[iPlace];  // multiply total by places's range
        nPacked += perm.b[iPlace];  // add place to total
    }
    return nPacked;
}

void CBalaGray::MakePerms(int nPlaces, const uint8_t *arrRange)
{
    m_nPlaces = nPlaces;
    m_arrRange.resize(nPlaces);
    int nPerms = 1;
    for (int iPlace = 0; iPlace < nPlaces; iPlace++) {  // for each place
        assert(arrRange[iPlace] > 1);   // radix must be at least binary
        m_arrRange[iPlace] = arrRange[iPlace];  // store range
        nPerms *= arrRange[iPlace]; // update range
    }
    m_arrPerm.resize(nPerms);
    for (int iPerm = 0; iPerm < nPerms; iPerm++) {
        PERM    perm;
        perm.dw = 0;
        int nVal = iPerm;
        for (int iPlace = 0; iPlace < nPlaces; iPlace++) {  // for each place
            int nRange = m_arrRange[iPlace];
            perm.b[iPlace] = nVal % nRange;
            nVal /= nRange;
        }
        m_arrPerm[iPerm] = perm;
    }
}

void CBalaGray::MakeGrayTable()
{
    // Build 2D table of permutations reachable from each permutation.
    // One row for each permutation, one column for each Gray neighbor.
    // Each element is a permutation index, and must be dereferenced.
    int nPlaces = m_nPlaces;
    int nGrayPerms = 0;
    for (int iPlace = 0; iPlace < nPlaces; iPlace++) {  // for each place
        nGrayPerms += m_arrRange[iPlace] - 1;   // one less than place's range
    }
    // Compute stride of Gray permutations array; to avoid multiplication,
    // round up stride to nearest power of two and convert it to a shift.
    unsigned long   iFirstBitPos;
    _BitScanReverse(&iFirstBitPos, nGrayPerms - 1);
    int nStrideShift = 1 << iFirstBitPos;
    m_arrGray.resize(m_arrPerm.size() << nStrideShift);
    int nPerms = GetPermCount();
    for (int iPerm = 0; iPerm < nPerms; iPerm++) {  // for each permutation
        int iCol = 0;
        PERM    rowPerm, colPerm;
        rowPerm.dw = m_arrPerm[iPerm].dw;
        for (int iPlace = 0; iPlace < nPlaces; iPlace++) {  // for each place
            int nRange = m_arrRange[iPlace];    // place's range
            for (int iVal = 0; iVal < nRange; iVal++) { // for of place's values
                if (iVal != rowPerm.b[iPlace]) {    // if value differs from row value
                    colPerm.dw = rowPerm.dw;    // column permutation is same as row
                    colPerm.b[iPlace] = iVal;   // except one place differs (Gray)
                    m_arrGray[(iPerm << nStrideShift) + iCol] = Pack(colPerm);
                    iCol++; // next column
                }
            }
        }
    }
    m_nGrayPerms = nGrayPerms;  // save in member var
    m_nGrayStrideShift = nStrideShift;
}

void CBalaGray::DumpPerm(const PERM& perm) const
{
    printf("[");
    for (int iPlace = 0; iPlace < m_nPlaces; iPlace++) {    // for each place
        printf("%d ", perm.b[iPlace]);
    }
    printf("]");
}

void CBalaGray::DumpPerms() const
{
    int nPerms = GetPermCount();
    for (int iPerm = 0; iPerm < nPerms; iPerm++) {
        for (int iPlace = 0; iPlace < m_nPlaces; iPlace++) {    // for each place
            printf("%d ", m_arrPerm[iPerm].b[iPlace]);
        }
        printf("\n");
    }
}

void CBalaGray::DumpGrayTablePerms() const
{
    int nPerms = GetPermCount();
    for (int iPerm = 0; iPerm < nPerms; iPerm++) {  // for each permutation
        DumpPerm(m_arrPerm[iPerm]);
        printf(": ");
        for (int iGray = 0; iGray < m_nGrayPerms; iGray++) {    // for each Gray neighbor
            int iPerm2 = m_arrGray[(iPerm << m_nGrayStrideShift) + iGray];
            DumpPerm(m_arrPerm[iPerm2]);
        }
        printf("\n");
    }
}

void CBalaGray::DumpSet() const
{
    printf("[");
    for (int iPlace = 0; iPlace < m_nPlaces; iPlace++) {    // for each place
        printf("%d", m_arrRange[iPlace]);
    }
    printf("]\n");
}

void CBalaGray::WriteBalanceToLog(int nImbalance, int nMaxTrans, int nMaxSpan)
{
    m_fOut << "balance = " << nImbalance << ", maxtrans = " << nMaxTrans << ", maxspan = " << nMaxSpan << '\n';
}

void CBalaGray::WriteSequenceToLog(int iDepth)
{
    int nPerms = GetPermCount();
    for (int iPlace = 0; iPlace < m_nPlaces; iPlace++) {    // for each place
        for (int iPerm = 0; iPerm < nPerms; iPerm++) {
            m_fOut << int(m_arrPerm[m_arrState[iPerm].iPerm].b[iPlace]) << ' ';
        }
        m_fOut << '\n';
    }
    m_fOut << '\n';
}

__forceinline bool CBalaGray::IsGray(PERM p1, PERM p2) const
{
    // Returns true if the given permutations differ by exactly one place.
    bool    bDiff = false;
    int nPlaces = m_nPlaces;
    for (int iPlace = 0; iPlace < nPlaces; iPlace++) {  // for each place
        if (p1.b[iPlace] != p2.b[iPlace]) { // if places differ
            if (!bDiff) {   // if first difference
                bDiff = true;   // set flag
            } else {    // not first difference
                return false;   // not Gray; early out
            }
        }
    }
    return bDiff;
}

void CBalaGray::Calc(int nPlaces, const uint8_t *arrRange)
{
    assert(nPlaces >= 0 && nPlaces <= MAX_PLACES);
    Reset();
    MakePerms(nPlaces, arrRange);
    MakeGrayTable();
//  DumpPerms();
//  DumpGrayTablePerms();
    int nPermGrays = m_nGrayPerms;
    int nGrayStrideShift = m_nGrayStrideShift;
    DumpSet();
    int nPerms = GetPermCount();
    printf("nPlaces=%d\n", nPlaces);
    printf("nPerms=%d\n", nPerms);
    int nBestImbalance = INT_MAX;
    int nBestMaxTrans = INT_MAX;
    int nBestMaxSpan = INT_MAX;
    m_arrState.resize(nPerms);
    uint64_t    nPasses = 0;
    uint64_t    nPermUsedMask[2] = {0}; // need 128 bits, as number of permutations may exceed 64
    int iDepth = 2; // first two levels are constant to save time; all sequences start with 0, 1
    m_arrState[1].iPerm = 1;
    m_arrState[1].nTrans.b[0] = 1;
    nPermUsedMask[0] = 0x3;
    int nStartDepth = iDepth;
    while (1) {
        nPasses++;
        int iPrevPerm = m_arrState[iDepth - 1].iPerm;
        int iGray = m_arrState[iDepth].iGray;
        int iPerm = m_arrGray[(iPrevPerm << nGrayStrideShift) + iGray]; // optimized 2D table addressing
        int iUsedMask = iPerm >= ULONGLONG_BITS;    // index selects one of two 64-bit masks
        uint64_t    nPermMask = 1ull << (iPerm & (ULONGLONG_BITS - 1));
        if (!(nPermUsedMask[iUsedMask] & nPermMask)) {  // if this permutation hasn't been used yet on this branch
            m_arrState[iDepth].iPerm = iPerm;   // save permutation index on stack
            int nMaxTrans;
            PERM    nTransCounts;
            int nImbalance = ComputeBalance(iDepth, nMaxTrans, nTransCounts);
            if (iDepth < nPerms - 1) {  // if incomplete sequence
#if DO_PRUNING
                // these constants may require tuning, see notes below
//              if (nMaxTrans > PRUNE_MAXTRANS || nImbalance > PRUNE_IMBALANCE) {   // slightly faster
                if (nImbalance > PRUNE_IMBALANCE) {
                    goto lblPrune;  // abandon this branch
                }
#endif
                // crawl one level deeper
                nPermUsedMask[iUsedMask] |= nPermMask;  // mark this permutation as used
                m_arrState[iDepth].nTrans.dw = nTransCounts.dw; // save current transition counts on stack
                iDepth++;   // increment depth to next permutation
                m_arrState[iDepth].iGray = 0;   // reset index of Gray neighbors
                m_arrState[iDepth].iPerm = 0;   // reset permutation index
                continue;   // equivalent to recursion, but less overhead
            } else {    // reached a leaf: complete sequence, a potential winner
                // if branch doesn't wrap around Gray
                if (!IsGray(m_arrPerm[m_arrState[0].iPerm], m_arrPerm[m_arrState[nPerms - 1].iPerm])) {
                    goto lblPrune;  // abandon this branch
                }
                // if max transition count or imbalance are worse than our current bests
                if (nMaxTrans > nBestMaxTrans || nImbalance > nBestImbalance) {
                    goto lblPrune;  // abandon this branch
                }
                int nMaxSpan = ComputeMaxSpan(iDepth);  // compute maximum span length
                // if max transition count and imbalance equal our current bests
                if (nMaxTrans == nBestMaxTrans && nImbalance == nBestImbalance) {
                    if (nMaxSpan >= nBestMaxSpan) { // if max span didn't improve
                        goto lblPrune;  // abandon this branch
                    }
                }
                // we have a winner, until something better comes along
                nBestMaxTrans = nMaxTrans;  // update best max transition count
                nBestImbalance = nImbalance;    // update best imbalance
                nBestMaxSpan = nMaxSpan;    // update best maximum span length
                printf("balance = %d, maxtrans = %d, maxspan = %d\n", nImbalance, nMaxTrans, nMaxSpan);
                WriteBalanceToLog(nImbalance, nMaxTrans, nMaxSpan);
                WriteSequenceToLog(iDepth);
            }
        }
        m_arrState[iDepth].iGray++; // increment Gray neighbor index
        if (m_arrState[iDepth].iGray >= nPermGrays) {   // if no more Gray neighbors for this permutation
lblPrune:
            if (iDepth <= nStartDepth) {    // if we're at same level where we started
                break;  // exit main loop
            } else {    // sufficient levels remain above us
                iDepth--;   // back up a level
                // restore bitmask that keeps track of which permutations we've used on this branch
                int iPerm = m_arrState[iDepth].iPerm;   // number of permutations may exceed 64
                int iUsedMask = iPerm >= ULONGLONG_BITS;    // index selects one of two 64-bit masks
                uint64_t    nPermMask = 1ull << (iPerm & (ULONGLONG_BITS - 1));
                nPermUsedMask[iUsedMask] &= ~nPermMask; // mark this permutation as available again
                m_arrState[iDepth].iGray++; // increment was skipped by continue statement above
                if (m_arrState[iDepth].iGray >= nPermGrays) {   // if no more Gray neighbors
                    goto lblPrune;  // keep backing up
                }
            }
        }
    }
    printf("done!\n");
}

__forceinline int CBalaGray::ComputeBalance(int iDepth, int& nMaxTrans, PERM& nTransCounts) const
{
    int nPlaces = m_nPlaces;
    PERM    nTrans;
    nTrans.dw = m_arrState[iDepth - 1].nTrans.dw;   // load latest transition counts from stack
    // compare current state to previous state
    PERM    sPrev, sCur;
    sPrev.dw = m_arrPerm[m_arrState[iDepth - 1].iPerm].dw;
    sCur.dw = m_arrPerm[m_arrState[iDepth].iPerm].dw;
    for (int iPlace = 0; iPlace < nPlaces; iPlace++) {  // for each place
        if (sCur.b[iPlace] != sPrev.b[iPlace]) {    // if place transitioned
            nTrans.b[iPlace]++; // increment place's transition count
        }
    }
    nTransCounts = nTrans;  // order matters; counts passed back to caller must exclude wraparound
    // account for wraparound; compare current state to initial state, which is assumed to be zero
    for (int iPlace = 0; iPlace < nPlaces; iPlace++) {  // for each place
        if (sCur.b[iPlace]) {   // if place transitioned
            nTrans.b[iPlace]++; // increment place's transition count
        }
    }
    // now that we have latest transition counts, compute their min and max
    int nMin = nTrans.b[0]; // initialize min and max to first transition count
    int nMax = nTrans.b[0];
    for (int iPlace = 1; iPlace < nPlaces; iPlace++) {  // for each transition count, excluding first
        int n = nTrans.b[iPlace];
        if (n < nMin)   // if less than min
            nMin = n;   // update min
        if (n > nMax)   // if greater than max
            nMax = n;   // udpate max
    }
    nMaxTrans = nMax;
    return nMax - nMin; // return difference
}

__forceinline int CBalaGray::ComputeMaxSpan(int iDepth) const
{
    int arrSpan[MAX_PLACES];
    int arrFirstSpan[MAX_PLACES];
    for (int iPlace = 0; iPlace < m_nPlaces; iPlace++) {    // for each place
        arrSpan[iPlace] = 1;    // initial span length is one
        arrFirstSpan[iPlace] = 0;   // first span length not set
    }
    int nMaxSpan = 1;
    PERM    sFirst, sPrev;
    sFirst.dw = m_arrPerm[m_arrState[0].iPerm].dw;  // store first state
    sPrev.dw = sFirst.dw;
    for (int iState = 1; iState <= iDepth; iState++) {  // for each state, excluding first
        PERM    s;
        s.dw = m_arrPerm[m_arrState[iState].iPerm].dw;  // compare this state to previous state
        for (int iPlace = 0; iPlace < m_nPlaces; iPlace++) {    // for each place
            if (s.b[iPlace] != sPrev.b[iPlace]) {   // if place transitioned
                if (arrSpan[iPlace] > nMaxSpan) // if span length exceeds max
                    nMaxSpan = arrSpan[iPlace]; // update max span length
                if (!arrFirstSpan[iPlace])  // if first span length hasn't been set
                    arrFirstSpan[iPlace] = arrSpan[iPlace]; // save first span length
                arrSpan[iPlace] = 1;    // reset span length
            } else {    // place didn't transition
                arrSpan[iPlace]++;  // increment span length
            }
        }
        sPrev = s;  // update previous state
    }
    // wrap around from last to first state
    for (int iPlace = 0; iPlace < m_nPlaces; iPlace++) {    // for each place
        if (sFirst.b[iPlace] != sPrev.b[iPlace]) {  // if place transitioned
            if (arrSpan[iPlace] > nMaxSpan) // if span length exceeds max
                nMaxSpan = arrSpan[iPlace]; // update max span length
        } else {    // place didn't transition
            arrSpan[iPlace] += arrFirstSpan[iPlace];    // compute wrapped span length
            if (arrSpan[iPlace] > nMaxSpan) // if span length exceeds max
                nMaxSpan = arrSpan[iPlace]; // update max span length
        }
    }
    return nMaxSpan;
}

void test()
{
// All cases want PRUNE_IMBALANCE = 3 unless specified otherwise below.
// Pruning greatly reduces runtime, but the results may not be optimal.
// Proven means exited normally with pruning disabled (DO_PRUNING = 0).
//
//  const uint8_t arrRange[] = {2, 10};     // proven
//  const uint8_t arrRange[] = {3, 9};
//  const uint8_t arrRange[] = {4, 8};
//  const uint8_t arrRange[] = {5, 7};
//  const uint8_t arrRange[] = {6, 6};
//  const uint8_t arrRange[] = {2, 9};      // proven
//  const uint8_t arrRange[] = {3, 8};
//  const uint8_t arrRange[] = {4, 7};
//  const uint8_t arrRange[] = {5, 6};
//  const uint8_t arrRange[] = {2, 8};      // proven
//  const uint8_t arrRange[] = {3, 7};
//  const uint8_t arrRange[] = {4, 6};
//  const uint8_t arrRange[] = {5, 5};
//  const uint8_t arrRange[] = {2, 7};      // proven
//  const uint8_t arrRange[] = {3, 6};      // proven
//  const uint8_t arrRange[] = {4, 5};      // proven
//  const uint8_t arrRange[] = {2, 6};      // proven
//  const uint8_t arrRange[] = {3, 5};      // proven
//  const uint8_t arrRange[] = {4, 4};      // proven
//  const uint8_t arrRange[] = {2, 5};      // proven
//  const uint8_t arrRange[] = {3, 4};      // proven
//  const uint8_t arrRange[] = {2, 4};      // proven
//  const uint8_t arrRange[] = {3, 3};      // proven
//  const uint8_t arrRange[] = {2, 3};      // proven
//  const uint8_t arrRange[] = {2, 2};      // proven
//  const uint8_t arrRange[] = {2, 2, 8};
//  const uint8_t arrRange[] = {2, 3, 7};
//  const uint8_t arrRange[] = {2, 4, 6};
//  const uint8_t arrRange[] = {2, 5, 5};
//  const uint8_t arrRange[] = {3, 3, 6};
//  const uint8_t arrRange[] = {3, 4, 5};
//  const uint8_t arrRange[] = {4, 4, 4};
//  const uint8_t arrRange[] = {2, 2, 7};
//  const uint8_t arrRange[] = {2, 3, 6};
//  const uint8_t arrRange[] = {2, 4, 5};
//  const uint8_t arrRange[] = {3, 3, 5};
//  const uint8_t arrRange[] = {3, 4, 4};
//  const uint8_t arrRange[] = {2, 2, 6};       // proven
//  const uint8_t arrRange[] = {2, 3, 5};
//  const uint8_t arrRange[] = {2, 4, 4};
//  const uint8_t arrRange[] = {3, 3, 4};
//  const uint8_t arrRange[] = {2, 2, 5};       // proven
//  const uint8_t arrRange[] = {2, 3, 4};       // proven
//  const uint8_t arrRange[] = {3, 3, 3};
//  const uint8_t arrRange[] = {2, 2, 4};       // proven
//  const uint8_t arrRange[] = {2, 3, 3};       // proven
//  const uint8_t arrRange[] = {2, 2, 3};       // proven
//  const uint8_t arrRange[] = {2, 2, 2};       // proven
//  const uint8_t arrRange[] = {2, 2, 2, 6};
//  const uint8_t arrRange[] = {2, 2, 3, 5};    // slow
//  const uint8_t arrRange[] = {2, 2, 4, 4};
//  const uint8_t arrRange[] = {2, 3, 3, 4};    // slow; wants PRUNE_IMBALANCE = 4
    const uint8_t arrRange[] = {3, 3, 3, 3};    // slow
//  const uint8_t arrRange[] = {2, 2, 2, 5};
//  const uint8_t arrRange[] = {2, 2, 3, 4};
//  const uint8_t arrRange[] = {2, 3, 3, 3};
//  const uint8_t arrRange[] = {2, 2, 2, 4};
//  const uint8_t arrRange[] = {2, 2, 3, 3};    // slow
//  const uint8_t arrRange[] = {2, 2, 2, 3};    // proven
//  const uint8_t arrRange[] = {2, 2, 2, 2};    // proven
//
// *** following cases require MORE_PLACES to be non-zero ***
//
//  const uint8_t arrRange[] = {2, 2, 2, 2, 4}; // wants PRUNE_IMBALANCE = 2
//  const uint8_t arrRange[] = {2, 2, 2, 3, 3}; // wants PRUNE_IMBALANCE = 4
//  const uint8_t arrRange[] = {2, 2, 2, 2, 3}; // wants PRUNE_IMBALANCE = 2
//  const uint8_t arrRange[] = {2, 2, 2, 2, 2};
//  const uint8_t arrRange[] = {2, 2, 2, 2, 2, 2};
//
    CBalaGray   bg;
    bg.Calc(_countof(arrRange), arrRange);
    fgetc(stdin);
}

int _tmain(int argc, _TCHAR* argv[])
{
    test();
    return 0;
}

Answer 1

Within CPLEX, I would use CPOptimizer.

For instance, in OPL

// 2 2 2 2



using CP;



int Size=4;
int r[1..Size]=[2,2,2,2];

int States=prod(i in 1..Size) r[i];

int fig[1..States][1..Size];

execute
{
  var index=0;
  for(var f1=1;f1<=r[1];f1++)
  for(var f2=1;f2<=r[2];f2++)
  for(var f3=1;f3<=r[3];f3++)
  for(var f4=1;f4<=r[4];f4++)
  {
    index++;
    fig[index][1]=f1;
    fig[index][2]=f2;
    fig[index][3]=f3;
    fig[index][4]=f4;
  }
}

dvar int x[1..States] in 1..States; // list of States in the right order
dvar int change[1..States] in 1..Size; // the figure that is different next time

dexpr int nbChanges[i in 1..Size]=count(change,i);
dexpr int inbalance=max(i in 1..Size) nbChanges[i]-min(i in 1..Size) nbChanges[i];

dvar int+ nochangeForThatManyTimes[1..States][1..Size] in 1..maxint;
dexpr int maxspan=max(i in 1..States,j in 1..Size)  nochangeForThatManyTimes[i][j];

minimize staticLex(inbalance,maxspan);
subject to
{
  x[1]==1;
  
  allDifferent(x);
  
  // Gray
  forall(i in 1..States,j in 1..Size) 
  ((fig[x[i]][j]==fig[x[(i<States)?(i+1):1]][j])==(j!=change[i]));
  
  
  forall(i in 2..States,j in 1..Size)  
  {
    (j==change[i-1]) =>   (nochangeForThatManyTimes[i][j]==1);
     (j!=change[i-1]) =>   (nochangeForThatManyTimes[i][j]==1+nochangeForThatManyTimes[i-1][j]);
  }



forall(j in 1..Size)  
(j==change[States]) ==
     (nochangeForThatManyTimes[1][j]==1);
 
   
}

execute
{
for(var i=1;i<=States;i++) 
{
  for(var j=1;j<=Size;j++) write(fig[x[i]][j]-1);
  writeln();
}  
writeln();
writeln("inbalance = ",inbalance);
writeln("maxspan = ",maxspan);
}

gives

0000
1000
1010
1011
1001
1101
0101
0001
0011
0010
0110
0111
1111
1110
1100
0100

inbalance = 0
maxspan = 6

and with 3,3,3,3 and a 12000 time limit I got

 OBJECTIVE: 1; 14
0000
1000
2000
2100
2110
2010
2012
0012
0010
0011
2011
1011
1010
1012
1022
1002
0002
0102
0202
0212
0210
2210
2200
2220
0220
0222
1222
2222
2202
2212
2112
1112
0112
0110
0120
0020
0021
0022
0122
0121
0111
0211
0221
1221
1220
1210
1211
1212
1202
1102
1122
1120
1110
1111
1121
1021
1020
2020
2120
2122
2102
2002
2022
2021
2121
2221
2211
2111
2101
2001
1001
0001
0101
1101
1201
2201
0201
0200
1200
1100
0100

inbalance = 1
maxspan = 14

and after 10 hours

OBJECTIVE: 1; 10
0000
0001
0002
0012
0010
0011
0111
0121
0221
1221
1021
0021
0022
0222
0122
0102
0112
0212
0202
0201
0101
1101
1001
1011
1012
1112
1110
0110
2110
2010
2020
2000
2002
1002
1022
1020
1120
0120
2120
2121
2221
2222
2202
2200
0200
1200
1000
1010
1210
2210
0210
0211
1211
1201
1202
1102
1100
0100
2100
2101
2102
2112
2012
2022
2122
1122
1121
1111
2111
2011
2021
2001
2201
2211
2212
1212
1222
1220
2220
0220
0020

inbalance = 1
maxspan = 10

In order to get

inbalance = 1
maxspan = 9

with

0000
1000
1001
1201
1200
1210
0210
2210
2211
2201
2200
2000
2010
1010
0010
0110
0112
0102
0202
0212
0012
1012
1011
1111
0111
0101
1101
1100
1102
1002
1022
1020
1120
1121
1021
0021
0011
0211
0221
0222
0122
1122
2122
2102
2202
1202
1222
1212
1211
1221
2221
2222
2212
2012
2022
0022
0002
2002
2001
2101
2121
0121
0120
0020
0220
1220
2220
2020
2021
2011
2111
2112
1112
1110
2110
2120
2100
0100
0200
0201
0001

I slightly improved the model

execute
{
  cp.param.timelimit=36000;
 
}

using CP;



int Size=4;
int r[1..Size]=[3,3,3,3];
int maxr=max(i in 1..Size) r[i];

int States=prod(i in 1..Size) r[i];

int fig[1..States][1..Size];

int which[i1 in 1..r[1]][i2 in 1..r[2]][i3 in 1..r[3]][i4 in 1..r[4]];

execute
{
  var index=0;
  for(var f1=1;f1<=r[1];f1++)
  for(var f2=1;f2<=r[2];f2++)
  for(var f3=1;f3<=r[3];f3++)
  for(var f4=1;f4<=r[4];f4++)
  {
    index++;
    fig[index][1]=f1;
    fig[index][2]=f2;
    fig[index][3]=f3;
    fig[index][4]=f4;
    which[f1][f2][f3][f4]=index;
  }
}



dvar int x[1..States] in 1..States; // list of States in the right order

dvar int y[1..States] in 1..States;


dvar int change[1..States] in 1..Size; // the figure that is different next time
dvar int move[1..States] in 1..maxr;

dexpr int nbChanges[i in 1..Size]=count(change,i);
dexpr int inbalance=max(i in 1..Size) nbChanges[i]-min(i in 1..Size) nbChanges[i];

dvar int+ nochangeForThatManyTimes[1..States][1..Size] in 1..maxint;
dexpr int maxspan=max(i in 1..States,j in 1..Size)  nochangeForThatManyTimes[i][j];

minimize staticLex(inbalance,maxspan);

subject to
{
  
//   inverse(x,y);
//  allDifferent(y);
  
  x[1]==1;
  forall(i in 1..States) move[i]<=r[change[i]];
  //inverse(x,y);
  
  change[1]==1;
  
  
  allDifferent(x);
  
  // Gray
//  forall(i in 1..States,j in 1..Size) 
//  {
//  (j!=change[i]) == (fig[x[i]][j]==fig[x[(i<States)?(i+1):1]][j]);
//  (j==change[i]) == ((fig[x[i]][j]+move[i]-1) mod r[j]+1==fig[x[(i<States)?(i+1):1]][j]);
//}  
 
 forall(i in 1..States)
   
   x[(i<States)?(i+1):1]
   ==which
   [(fig[x[i]][1]+(1==change[i])*move[i]-1) mod r[1]+1]
   [(fig[x[i]][2]+(2==change[i])*move[i]-1) mod r[2]+1]
   [(fig[x[i]][3]+(3==change[i])*move[i]-1) mod r[3]+1]
   [(fig[x[i]][4]+(4==change[i])*move[i]-1) mod r[4]+1]
   ;
   
   
 
  
  
  inferred(change);
  inferred(move);
  inferred(nochangeForThatManyTimes);
 
 
  inbalance>=States mod 2;
  
  forall(i in 2..States,j in 1..Size)  
  {
    (j==change[i-1]) =>   (nochangeForThatManyTimes[i][j]==1);
     (j!=change[i-1]) =>   (nochangeForThatManyTimes[i][j]==1+nochangeForThatManyTimes[i-1][j]);
  }



forall(j in 1..Size)  
(j==change[States]) ==
     (nochangeForThatManyTimes[1][j]==1);
 
   
}

execute
{
for(var i=1;i<=States;i++) 
{
  for(var j=1;j<=Size;j++) write(fig[x[i]][j]-1);
  writeln();
}  
writeln();
writeln("inbalance = ",inbalance);
writeln("maxspan = ",maxspan);
}

NB: You can use CPLEX for free in the cloud with this OPL API

Answer 2

Well, after some tinkering, I have an Integer Program running for this, that I think is producing quality results. Tried a couple approaches...each had differing limitations

It is a little grotesque in parts as the counting of repeat digits is quite cumbersome.

It really bogs down for things with ~30 states or more, so it's not going to make it to the finish line. :) I think it is much more nimble if I remove the repeat counting, and I'll tinker a bit more. In the interim, here are some results for the cases not marked as proven on your web page. The (4, 6) run (second run) is an improvement, the other 2 are now "proven" as stated, perhaps with a different sequence, I didn't x-check.

I'll update later with any other improvements.

starting run:  (3, 7)
WARNING: Initializing ordered Set PR_flat with a fundamentally unordered data
    source (type: set).  This WILL potentially lead to nondeterministic
    behavior in Pyomo

Problem: 
- Name: unknown
  Lower bound: 1.3
  Upper bound: 1.3
  Number of objectives: 1
  Number of constraints: 3744
  Number of variables: 3539
  Number of binary variables: 3557
  Number of integer variables: 3560
  Number of nonzeros: 3
  Sense: minimize
Solver: 
- Status: ok
  User time: -1.0
  System time: 353.01
  Wallclock time: 305.85
  Termination condition: optimal
  Termination message: Model was solved to optimality (subject to tolerances), and an optimal solution is available.
  Statistics: 
    Branch and bound: 
      Number of bounded subproblems: 1018
      Number of created subproblems: 1018
    Black box: 
      Number of iterations: 667719
  Error rc: 0
  Time: 306.00031781196594
Solution: 
- number of solutions: 0
  number of solutions displayed: 0

11
01
02
00
10
20
21
22
12
13
23
03
05
25
26
24
14
04
06
16
15
max imbalance: 1.0
max repeats: 3.0
starting run:  (4, 6)
WARNING: Initializing ordered Set PR_flat with a fundamentally unordered data
    source (type: set).  This WILL potentially lead to nondeterministic
    behavior in Pyomo

Problem: 
- Name: unknown
  Lower bound: 0.2
  Upper bound: 0.2
  Number of objectives: 1
  Number of constraints: 4854
  Number of variables: 4619
  Number of binary variables: 4640
  Number of integer variables: 4643
  Number of nonzeros: 3
  Sense: minimize
Solver: 
- Status: ok
  User time: -1.0
  System time: 34.21
  Wallclock time: 34.89
  Termination condition: optimal
  Termination message: Model was solved to optimality (subject to tolerances), and an optimal solution is available.
  Statistics: 
    Branch and bound: 
      Number of bounded subproblems: 1
      Number of created subproblems: 1
    Black box: 
      Number of iterations: 14167
  Error rc: 0
  Time: 34.923232078552246
Solution: 
- number of solutions: 0
  number of solutions displayed: 0

10
13
33
34
14
15
35
32
02
03
23
22
12
11
01
00
30
31
21
24
04
05
25
20
max imbalance: 0.0
max repeats: 2.0
starting run:  (5, 5)
WARNING: Initializing ordered Set PR_flat with a fundamentally unordered data
    source (type: set).  This WILL potentially lead to nondeterministic
    behavior in Pyomo

Problem: 
- Name: unknown
  Lower bound: 1.3
  Upper bound: 1.3
  Number of objectives: 1
  Number of constraints: 5256
  Number of variables: 5011
  Number of binary variables: 5033
  Number of integer variables: 5036
  Number of nonzeros: 3
  Sense: minimize
Solver: 
- Status: ok
  User time: -1.0
  System time: 915.71
  Wallclock time: 634.99
  Termination condition: optimal
  Termination message: Model was solved to optimality (subject to tolerances), and an optimal solution is available.
  Statistics: 
    Branch and bound: 
      Number of bounded subproblems: 1764
      Number of created subproblems: 1764
    Black box: 
      Number of iterations: 1855323
  Error rc: 0
  Time: 635.0473001003265
Solution: 
- number of solutions: 0
  number of solutions displayed: 0

11
01
31
33
34
44
04
03
00
40
41
42
22
23
43
13
12
02
32
30
10
20
21
24
14
max imbalance: 1.0
max repeats: 3.0