Proof involving the Birkhoff polytope Setup:
Fix $n$. Let $\mathcal{B}$ denote the $n$-dimensional Birkhoff polytope, i.e. the set of all $n \times n$ doubly stochastic matrices. There is a fixed natural number $N(n)$. We drop the argument $n$ of $N$ for the rest of the question.
Suppose $\Gamma=\mathcal{B}^N$, i.e. the set of all ordered sets of $N$ points, each chosen independently from the Birkhoff polytope. Let a typical element of $\Gamma$ be denoted by $x \equiv (x^1,\cdots,x^N)$, where $x^1,\cdots,x^N$ are bistochastic matrices. Let $x^k_{ij}$ denote the $(i,j)$-th entry of the bistochastic matrix $x^k$.
Now we define $S=\{x \in \Gamma: Ax=0, Bx \geq 0\}$, where $A$ and $B$ are given linear transformations ($A$ and $B$ depend on $n$). Each row of each of $A$ and $B$ is of the form $[0, \cdots, 0, 1, -1, 0, \cdots 0]$, i.e. each linear restriction imposed by $A$ and $B$ only involve difference between two terms. It may also be worth mentioning that each of the difference-based restrictions contained in $A$ ($B$) is of the form $x^k_{ij}-x^l_{ij}=0\;(\geq 0)$, i.e. difference between corresponding entries of two bistochastic matrices from $\Gamma$. Other than these restrictions, $A$ and $B$ can be completely general.
What I'm trying to prove:
$S$ is obviously a polytope. I'm trying to prove it has integral extreme points, i.e. its extreme points consist only of entries in $\{0,1\}$. It is my conjecture that this holds for any $A$ and $B$ of the form described above.
My approach: 
Consider $I=\{x \in \Gamma: Ax=0\}$. Clearly $S \subseteq I$. Then I have used the exact same logic as given in the proof of the Birkhoff-von Neumann theorem here to argue that $I$ cannot have fractional extreme points. Essentially the logic is as follows: If there exists any entry of any element of $\Gamma$ which is in $(0,1)$, then there has to exist a cycle - of even length - of entries in $(0,1)$, such that each of the entries in this cycle can be increased and decreased by the same $\epsilon$ amount, while still staying inside $\Gamma$. In my case I know the cycle has to be of even length because the part of the cycle inside each bistochastic matrix has to be of even length, by the Birkhoff-von Neumann proof, and each point $x \in \Gamma$ is just an ordered set of bistochastic matrices. Therefore $I$ has integral extreme points.
By the above logic, any subset of $\Gamma$, defined by equalities of the form $x^k_{ij}-x^l_{ij}=0$ has extreme points consisting only of entries in $\{0,1\}$. 
Questions:
First question - is this argument correct?
Secondly, how do we bring the inequality constraints, $Bx \geq 0$, into this framework and make any conclusions about integral or fractional extreme points?
Thank you for your help.
EDIT1: 
Subsequently I've thought as follows (basically continuing the above approach but trying to incorporate the $Bx \geq 0$ constraints). Suppose there exists a point $x \in S$ with a fractional element somewhere. Suppose at this point some of the $Bx \geq 0$ inequalities hold with equality (zero of the $Bx \geq 0$ inequalities can also hold with equality). Then just include those equalities in the set of equality constraints in $I$ and repeat the same argument. The $Bx \geq 0$ inequalities which do not bind don't pose any difficulty because by finiteness of the system (number of coordinates and number of constraints) we can always choose $\epsilon$ small enough to make the proof work. 
Again, I'm not sure if this is correct. Thanks in advance everyone.  
EDIT2:
If you want to impose additional assumptions on $A$ or $B$ (but they have to satisfy the properties already specified) that's fine too, for the time being. 
 A: I cannot follow your logic. Why do you get a full cycle? It may well happen that some entry in a matrix $x^i$ can be increased/decreased, but no other entries in its row or column have this freedom, since they are constrained by conditions coming from $A$ or $B$ that link them to corresponding entries of other matrices $x^j$.
I also think the claim of your question is false, even if $B$ has no rows (so the $n$-tuple is only restricted by equalities between some entries, not by inequalities). Here is a counterexample: Set $n = 3$ and $N = 2$, and let our two bistochastic matrices be $x^1 = \begin{pmatrix} a & b & c \\ d & e & f \\ g & h & i \end{pmatrix}$ and $x^2 = \begin{pmatrix} A & B & C \\ D & E & F \\ G & H & I \end{pmatrix}$ (sorry for re-using the letters $A$ and $B$). Further constrain them by requiring $f = F$ and $g = G$. (This yields a $2$-row matrix $A$ and a $0$-row matrix $B$.) Then, the resulting polytope (in an $18$-dimensional space) has $(1/2, 0, 1/2, 0, 1/2, 1/2, 1/2, 1/2, 0, 1/2, 1/2, 0, 0, 1/2, 1/2, 1/2, 0, 1/2)$ as one of its vertices. This vertex corresponds to the two matrices $x^1 = \begin{pmatrix} 1/2 & 0 & 1/2 \\ 0 & 1/2 & 1/2 \\ 1/2 & 1/2 & 0 \end{pmatrix}$ and $x^2 = \begin{pmatrix} 1/2 & 1/2 & 0 \\ 0 & 1/2 & 1/2 \\ 1/2 & 0 & 1/2 \end{pmatrix}$. There is a Sudoku-style manual proof of the fact that this is indeed a vertex; alternatively, here is a SageMath program that computes all the $18$ vertices:
Q.<a,b,c,d,e,f,g,h,i,A,B,C,D,E,F,G,H,I> = PolynomialRing(QQ)

vars = Q.gens()

def ieq(p):
    # Transform a degree-1 polynomial like `a + 2b - c - 5`
    # into a tuple of coefficients that can later be
    # entered into the ``ieqs`` or ``eqns`` attributes of
    # the ``Polyhedron`` constructor, where it will be
    # interpreted as `a + 2b - c \geq 5` or
    # `a + 2b - c = 5`, respectively.
    return map(QQ, tuple([p.coefficient({j: 0 for j in vars})] + [p.coefficient({j: int(j==i) for j in vars}) for i in vars]))

birk1ieqs = [a, b, c, d, e, f, g, h, i]
birk1eqs = [a+b+c-1, d+e+f-1, g+h+i-1, a+d+g-1, b+e+h-1, c+f+i-1]
birk2ieqs = [A, B, C, D, E, F, G, H, I]
birk2eqs = [A+B+C-1, D+E+F-1, G+H+I-1, A+D+G-1, B+E+H-1, C+F+I-1]
extraeqs = [f-F, g-G]

full_ieqs = map(ieq, birk1ieqs) + map(ieq,birk2ieqs)
full_eqs = map(ieq, birk1eqs) + map(ieq,birk2eqs) + map(ieq,extraeqs)

print Polyhedron(ieqs=full_ieqs, eqns=full_eqs).vertices()

