Espero te ayude

from itertools import product
from math import ceil
numeros = [0, 1, 6, 8, 9]
nro_mezcla=[]
n=3 # puedes usar cualquier valor de n 
for i in product(numeros,repeat=n):
    aux =ceil(n/2)
    c=0
    for j in range(0,aux):
        # aprovechando que puedo avanzar positivamente y negativamente 
        if i[-(j+1)]+i[j]==15: # si la suma de 6 y 9 es 15
            c+=1
        elif i[-(j+1)]+i[j]==16: # si la suma de 8 y 8 es 16
            c+=1
        elif i[-(j+1)]+i[j]==0: # si la suma de 0 y 0 es 0 XD 
            c+=1
        elif i[-(j+1)]+i[j]==2: # Si la suma de 1 y 1 es 2
            c+=1
    if c==aux:
        nro_mezcla.append(i)
print(nro_mezcla)

# Disfruta programando 

Mi solución usa un generador recursivo.

La idea es poner en el centro un valor de partida. Puede ser un único dígito (si n es impart) o un par estreboscopico de dos elementos. Luego se hace crecer el valor agregando por la izquierda y derecha las mitades de cada par estroboscopico elemental.

digits = ['0', '1', '8']
pairs = ['00', '11', '69', '88', '96']

def estro(n):

    if n % 2:
        medio = n // 2
        for item in estro(n-1):            
            left = item[:medio]
            right = item[medio:]
            for d in digits:
                yield "".join([left, d, right])            
    elif n == 2:
        for p in pairs:
            yield p
    elif n > 2:
        for item in estro(n - 2):
            for p in pairs:     
                yield "".join([p[0], item, p[1]])
    else:
        yield ""

La ùnica complicación es que también devuelve los valores que comienzan con cero, por lo que necesito otra función para filtrarlos:

def get_estro(n):
    for value in estro(n):
        if value[0] != '0':
            yield value

y para probar

for x in get_estro(4):
    print(x, end=", ")
print()

produce:

candid@dell ~ $ python3 estro.py
1001, 6009, 8008, 9006, 1111, 6119, 8118, 9116, 1691, 6699, 8698, 9696, 1881, 6889, 8888, 9886, 1961, 6969, 8968, 9966, 

Lo que necesitas para generar los números no son permutaciones, es el producto cartesiano de los valores de la lista:

>>> list(itertools.permutations([1, 2, 3], 2))
[(1, 2), (1, 3), (2, 1), (2, 3), (3, 1), (3, 2)]
>>> list(itertools.product([1, 2, 3], repeat=2))
[(1, 1), (1, 2), (1, 3), (2, 1), (2, 2), (2, 3), (3, 1), (3, 2), (3, 3)]

No se si este será el código más eficiente para calcular esos números, pero lo que se me ocurre es que se recorra cada elemento del iterable resultante. Si el primer elemento es 0, se omite (continue) Si la cantidad de elementos es impar, entonces el elemento central debe ser 0, 1, 8. Si no se cumple, se omite (continue)

Ahora en el for se evalua si el elemento[i] corresponde a su 'opuesto' elemento[longitud-i]

(1, 6, 8, 9, 1)
 ^  ^  *  ^  ^
 |  +-----+  |
 +-----------+

Y finalmente la tupla que cumple con los requisitos se añade a la lista arreglop. El uso de map(str, tup) es porque join() no acepta valores int.

import itertools

numeros = [0, 1, 6, 8, 9]
n = 3
prod = itertools.product(numeros, repeat=n)

arreglop = []
for tup in prod:
    if tup[0] == 0:
        continue
    if (n % 2 > 0) and (tup[int(n/2)] not in [0, 1, 8]):
        continue
    
    estrob = True
    for i in range(int(n/2)):
        if (tup[i] == 0) and (tup[n-i-1] == 0):
            pass
        elif (tup[i] == 1) and (tup[n-i-1] == 1):
            pass
        elif (tup[i] == 6) and (tup[n-i-1] == 9):
            pass
        elif (tup[i] == 8) and (tup[n-i-1] == 8):
            pass
        elif (tup[i] == 9) and (tup[n-i-1] == 6):
            pass
        else:
            estrob = False
            break
    if estrob:
        arreglop.append(''.join(map(str, tup)))

print(arreglop)

Se me ocurre una idea mejor que la de generar todas las combinaciones posibles para luego quedarse sólo con algunas, y sería generar directamente sólo aquellas combinaciones que sabemos que son "estroboscópicas". ¿Cómo? Ahora explico la idea.

Si n es par

Si el número de cifras es par, la cosa es relativamente sencilla. Por ejemplo, considera el caso n=2. Lo primero sería tomar la mitad de n, en este caso 1, y generar todas las combinaciones (incluso con repetición, por lo que no serían realmente las permutaciones) de las cifras 0,1,6,8,9 tomadas de n/2 en n/2, para seguidamente concatenar a cada una de esas combinaciones la combinación necesaria para que el número sea estroboscópico, pero saltándose aquellas que comiencen por cero.

Por ejemplo, la primera "combinación" con n=1 sería "0". Ya que comienza por "0" nos la saltamos. La siguiente sería "1". Le concatenamos al otro lado "1" que es la versión "girada" del 1. La siguiente sería "6". Le concatenamos al otro lado "9" que es su girado. La siguiente sería "8". Le concatenamos al otro lado "8", que sería su versión girada. Finalmente la última sería "9". Le concatenamos al otro lado "6".

Así obtenemos: "11", "69", "88", y "96", que son todas las que buscábamos.

Otro ejemplo para n=4. La mitad sería n/2=2. Las primeras combinaciones serían:

• "00". Empieza por 0, me la salto

• "01". Empieza por 0 me la salto

• "06". Me la salto,

  etc.. así hasta la primera que no empiece por 0

• "10". Le concateno "01" que es su versión girada y obtengo "1001"

• "11". Le concateno "11" y obtengo "1111".

• etc.

Bueno, creo que la idea ya se pilla. ¿Cómo implementarla? Necesitamos una función que nos vaya dando las combinaciones de esas cifras (itertools.product vale) y otra que nos dé la versión "girada" de una secuencia dada. Esta otra función sería tal que si le pasas la secuencia ("1", "0") te devuelve ("0", "1"). No es tan fácil como recorrer la tupla al revés, ya que si le pasas ("1", "6") te debe devolver ("9", "1").

La forma más simple de implementar esta segunda función sería tener un mapa (diccionario) de cada cifra admisible y cuál es su versión girada. Por tanto podría ser así:

mapa = {
    "0": "0",
    "1": "1",
    "6": "9",
    "8": "8",
    "9": "6"
}

def version_girada(numero):
  return tuple(mapa[c] for c in numero[::-1])

Observa cómo en la expresión generadora recorro la tupla numero "marcha atrás" con numero[::-1] y para cada cifra que obtengo, uso su versión girada según el mapa.

Haciendo uso de esta función, es posible implementar la idea que expuse más arriba para crear los números estroboscópicos en el caso de que n sea par.

from itertools import product

def numeros_estroboscopicos_n_par(n):
  resultado = []                      # Aqui iremos guardando los que salen
  mitad = n//2
  for parte in product(mapa, repeat=mitad): # Generar combinaciones de las cifras válidas
    if parte[0] == "0":                   # Nos saltamos los que empiezan por cero
      continue
    numero = parte + version_girada(parte)
    resultado.append("".join(numero))
  return sorted(resultado)

Podemos probar que funciona correctamente para n par:

>>> numeros_estroboscopicos_n_par(2)
['11', '69', '88', '96']

>>> numeros_estroboscopicos_n_par(4)
['1001', '1111', '1691', '1881', '1961', '6009', '6119', '6699', '6889',
 '6969', '8008', '8118', '8698', '8888', '8968', '9006', '9116', '9696',
 '9886', '9966']

Si n es impar

Por no darte el ejercicio completamente resuelto, te dejo que este los implementes tú, pero te doy la idea clave:

Si n es impar podemos comenzar generando todos los números estroboscópicos de longitud n-1. Y después, a cada uno de ellos le metemos en el centro otra cifra que sea simétrica.

Por ejemplo, para n=3. Hemos visto que para n=2 la lista salía "11", "69", "88", "96". Pues ahora a cada uno de ellos:

• Le metemos en el centro un 0 y salen: "101", "609", "808", "906"
• Le metemos en el centro un 1 y salen: "111", "619", "818", "916"
• Le metemos en el centro un 8 y salen: "181", "689", "888", "986"

Y de este modo ya las hemos generado todas.

Para implementar esto necesitas otra lista que contenga sólo los números simétricos, es decir "0", "1", "8", y el pseudocódigo sería entonces:

• Inicializar resultados con lista vacía
• Para cada número estroboscópico de longitud n-2
  • Para cada cifra simétrica
    • Añadir cifra en el centro del número estroboscópico
    • Añadir número resultante a resultados

Lo de añadir la cifra en el centro se puede hacer con slices: numero[:mitad] + cifra + numero[mitad:]

Juntando todo (cualquier n)

• Si n es impar, hacer m = n-1, si n es par hacer m=n
• Generar todos los numeros estroboscópicos de longitud m (par) con la función numeros_estroboscopicos_n_par(m).
• Si n era par, retornar el resultado
• En caso contrario usar la lista que se ha encontrado para insertar en el centro las cifras simétricas (como se explicó en el segundo apartado) y retornar esta otra lista.