An interesting problem. I've done a bit of work in texture compression and this sounds something like a generalisation of Campbell et al's "Color Cell Compression". It's also a little like a feature we were asked to include in the Dreamcast VQ compressor so that sub-palettes could be swapped to create different colour schemes on textures.
I was thus thinking that maybe you could adapt a standard VQ algorithm, applied in stages, to the problem.
You say that
"the exact colors, the arrangement of colors into subpalettes, and the subpalette for each tile of the image are all unknowns"
but can we assume the X*Y tile size is predetermined?
Assuming that is the case, perhaps the following might be worth trying:
1) Guess an initial palette size for the entire image, perhaps,
$$Psize=(M * SizeOfEachSubpalette)$$
where $1 < M < N$.
Use standard VQ quantisation to generate this palette and map the pixels to it.
2) Now consider each $X*Y$ tile to be a "super pixel" with $Psize$ channels. Set each channel value in the tile's "super pixel" to the number of pixels that mapped to the corresponding colour in our initial palette.
For example, if in our initial palette, entry $j$ is [0x11,0x21,0x31] and 5 pixels in the tile mapped to that colour, set the $j'th$ component in the tile's superpixel to 5.
Run a second VQ operation on the set of all "super pixels" using '$Psize$-dimensional space' to generate $N$ "palette" entries. (Actually, you don't need the 'palette' values - the VQ process is just meant as a means of partitioning the tiles into $N$ sets)
Hopefully, this will now identify which tiles should share the same sub-palette.
3) For each tile set, collect all the pixels in the tiles, and apply a 3rd VQ operation to identify the $SizeOfEachSubpalette$ colours needed for all tiles in that set.
I would think some sort of clever GLA/k-means process might be needed for step 2 to juggle the assignments of tiles to sub-palettes. The choice of M might also be interesting.
