Enquête Star avec optimisation des combinaisons

Qu'est-ce que c'est

• Enquêtez sur 1000 étoiles d'affilée avec une fusée à photons.
• Voyons s'il y a des extraterrestres. Cependant, toutes les étoiles ne peuvent pas être examinées.
• La durée de l'enquête varie d'une étoile à l'autre. La durée totale de l'enquête sur les étoiles à étudier doit être inférieure à 10 000 jours.
• La probabilité de découverte est estimée pour chaque étoile. Maximisez la somme des probabilités de découverte des étoiles étudiées.

Essayez de résoudre

Il peut être considéré comme un problème de sac à dos. Résolvons-le avec Python. Voir le lien de référence pour l'optimisation mathématique avec Python.

python3

import numpy as np
from pulp import *
np.random.seed(1)
Nombre d'étoiles= 1000
Heure de l'enquête= np.random.randint(10,100,Nombre d'étoiles)
Probabilité de découverte= np.random.random(Nombre d'étoiles)/100000
m = LpProblem(sense=LpMaximize)
x = [LpVariable('x%.4d'%i, cat=LpBinary) for i in range(Nombre d'étoiles)]
m += lpDot(Probabilité de découverte,x)
m += lpDot(Heure de l'enquête,x) <= 10000
m.solve()
print(value(m.objective)) #Somme des probabilités de découverte
>>>
0.0022822674119170536

En fait, les fusées ont une portée de croisière maximale. Ici, le nombre maximum de sauts + 1 est mx, pas simplement la distance. Voyons ce qui arrive à la fonction objectif lorsque nous changeons mx. L'axe horizontal est mx et l'axe vertical est la fonction objectif.

python3

r = []
for mx in range(4,17):
    m = LpProblem(sense=LpMaximize)
    x = [LpVariable('x%.4d'%i, cat=LpBinary) for i in range(Nombre d'étoiles)]
    m += lpDot(Probabilité de découverte,x)
    m += lpDot(Heure de l'enquête,x) <= 10000
    for i in range(Nombre d'étoiles-mx+1):
        m += lpSum(x[i:i+mx]) >= 1 #Rechercher au moins un emplacement dans mx
    m.solve()
    r.append(value(m.objective))

%matplotlib inline
import matplotlib.pyplot as plt
plt.plot(range(4,17),r)
plt.hlines(0.0022822674119170536,4,16);

Essayez d'agrandir.

python3

plt.plot(range(4,17),r)
plt.hlines(0.0022822674119170536,4,16)
plt.xlim((9,16))
plt.ylim((0.00227,0.0023));

Avec le solveur gratuit par défaut CBC, en raison d'erreurs, plus la solution est restrictive, mieux c'est.

Le solveur commercial a été résolu plus précisément.

Lien de référence

• Python dans l'optimisation --Qiita

c'est tout