Salut !

J'ai besoin d'un coup de main pour un problème d'optimisation. J'ai étudié le simplexe et la programmation linéaire il y a quelques années, mais j'ai l'impression que mon problème n'est pas résoluble de cette manière "directement" et qu'il faut le modifier un peu.

-
Edité par Transistored 25 juillet 2020 à 17:58:32

Merci de ta réponse,

j'insiste vraiment sur le fait de pouvoir généraliser "facilement", ce que j'ai du oublier d'écrire dans mon post initial. Et par généraliser j'entend pouvoir rajouter des recettes, passer d'un système avec 6 composants à plus, passer d'une ressource initiale à plusieurs etc.

Premièrement la fonction de scoring est assez arbitraire (et ne change pas fondamentalement le problème je pense).

Mais pour l'exemple disons qu'on veut : max -10P + 3A + 2R + F

j'ai implémenté un parcours d'arbre en largeur que je colle ci-dessous

initial=[480,0,0,0,0,0]
#P,A,H,L,F,R
a0=np.array([-30,20,10,0,0,0])
a1=np.array([0,20,0,-60,0,0])
a2=np.array([-30,0,20,0,0,20])
a3=np.array([-60,0,20,130,0,0])
a4=np.array([0,0,0,-40,0,20])
a5=np.array([0,-30,0,0,-30,60])
a6=np.array([-60,0,0,30,40,0])
a7=np.array([0,0,-60,0,40,0])
a=[a0,a1,a2,a3,a4,a5,a6,a7]

def possib(x,a):
    """
    teste s'il est possible de faire l'action a en étant dans l'état x
    """
    return int(sum(n < 0 for n in x+a)==0)
def score(x):
    return(##############)

def parcours(hist,etat):
    done=1
    nxt=[]
    #print(hist)
    for i in range(0,8):
        if possib(etat,a[i]):
            nxt.append(parcours(hist+[i],etat+a[i]))
            done=0
    if done:
        return [hist,etat]
    else:
        sc=[score(nxt[i][1]) for i in range(0,len(nxt))]       
        ind=np.argmax(sc)
        return(nxt[ind])

Pour 120P initiaux, on obtient 847 noeuds visités, pour 150P initiaux on est à presque 20000 noeuds, alors imagine avec 480...

C'est ce pourquoi l'idée du parcours d'arbre me gêne un peu. 

Pour ce qui est du temps de calcul grossier, j'ai considéré que dans le pire cas on avait une branche de longueur 16 (480/30 dans le cas ou on ferait 16 fois la recette 1). Chaque noeud a 8 fils dans le pire cas, donc 8^16. En prenant un processeur à 4Ghz et en supposant que le parcours d'un noeud soit en O(1) on aura environ 8^16/4e6 secondes à attendre, ce qui est assez long haha. Je sais que c'est loin d'être un bon calcul, mais l'approximation est assez bonne avec des valeurs de P entre 60 et 150.

Le problème avc ta démarche, c'est que tu prends en compte l'ordre.

Si on applique la recette ZA puis la recette ZC, ou si on applique la recette ZC, puis la recette ZA, tu comptes ça comme 2 scénarios différents. Alors qu'en fait, dans les 2 cas, on a transformé 60P en 20A+20R+30H.

Ce que tu peux faire, c'est à chaque fois tester si tu retrouves une configuration déjà analysée. 

Après ZA+ZC,  tu as en rayon 420P+20A+20R+30H

Quand tu envisages ZC+ZA, tu constates que tu tombes sur 420P+20A+20R+30H  ... et que ceci a déjà été analysé. Donc tu peux couper cette branche.

Si on n'avait que les 4 recettes qui partent de l'ingrédient P, on aurait environ une centaines de cas à analyser (cf mon 1er message), alors que toi, tu analyses environ 4^10 = 1 Million de cas.

Mais en terme de performance, tu gagneras encore si tu réussis à implémenter ce que je proposais.

Effectivement tu as raison, il n'est pas utile de s'occuper de l'ordre si on consomme la même ressource.

j'ai implémenté en ne prenant pas en compte l'ordre : on stocke non pas l'historique des recettes utilisées mais une liste Q telle que Q[i] te donne le nombre de fois ou tu as utilisé la recette i pour arriver à cet état.

Lorsqu'on visite un noeud, on regarde si cette matrice d'historique a déja été visitée. A mon avis cela revient à ce que tu disais, et ca ne prend pas en compte l'ordre. 

Malheuresement, bien qu'on gagne quand meme pas mal de temps pour des petites valeurs initiales de P, ca reste trop long pour Pi=480 (>10 minutes). Pour la curiosité je vais comparer le nombre de noeuds explorés dans les deux méthodes !

Je suis plus à la recherche d'une solution mathématique/optimisation que algorithmique en réalité ^^ 

Ici, sur cet exemple, voici un raisonnement qui a priori donne la meilleure solution.

La fonction d'évaluation est : maximiser -10P+3A+2R+F

Donc on peut dire P=-10, A=3, R=2 et F=1 ; en fait, ce sont des minorants. S'il s'avère qu'il y a une instruction 30F=2R, alors F vaut quelque chose entre 1 et 4/3

Ici, on a une instruction 40L=20R ; donc L=2R et L=1  ( sauf si malencontreusement, on avait seulement 10L par exemple)

L'instruction 60L=20A nous dit également L=1

Les 4 instructions à partir de P donnent :

30P=20A+10H=60+10H

30P=20R+20H=40+20H

60P=130L+20H=130+20H

60P=40F+30L=70

H vaut combien ? 60H=40F , donc H=2/3 si on a soit 60H, soit 120H, soit ...

Donc sur cet exemple, la meilleure solution, c'est de convertir nos 480H par l'instruction 60P=130L+20H, ce qui donne 1040L+160H

Les 1040L sont convertis via l'instruction 40L=20R

Et les 160H sont convertis via 60H=40F, donc

480P=1040L+160H=520R+80F+40H

Et on a un note de 520*2+80 = 1120 , et il nous reste nos 40H.

Comment généraliser le traitement au cas général ?

j'ai réussi à trouver une solution pour convertir ce problème en programme linéaire !

Supposons qu'on ait encore ces recettes, et P=480 initialement

(1) 3P -> 2 A + 1H

(5) 6L -> 2 A

Vient maintenant le temps de la fonction à maximiser !

Supposons qu'on dispose d'un score nA pour A, nF pour F et nR pour R

il s'agit tout simplement de max : nA*(2*x1+2*x5-3*x7) + nF*(4*x4-3*x7+4*x8) + nR*(2*x2+2*x3+2*x5+6*x7)

On obtient bien les mêmes valeurs qu'avec le parcours d'arbre, mais en 40ms. Voici le code pour cet exemple. J'ai pu généraliser assez facilement (pour pas avoir à calculer les équations à la main comme j'ai fais pour cet exemple) en utilisant le même principe. Je le mettrai probablement dans un github quand ca sera plus propre !

Merci de ton aide en tt cas :))

from __future__ import print_function
from ortools.linear_solver import pywraplp


def main():
    # Create the mip solver with the CBC backend.
    solver = pywraplp.Solver('simple_mip_program',
                         pywraplp.Solver.CBC_MIXED_INTEGER_PROGRAMMING)

    infinity = solver.infinity()
    valmax=20000
    x1 = solver.IntVar(0, valmax, 'x1')
    x2 = solver.IntVar(0, valmax, 'x2')
    x3 = solver.IntVar(0, valmax, 'x3')
    x4 = solver.IntVar(0, valmax, 'x4')
    x5 = solver.IntVar(0, valmax, 'x5')
    x6 = solver.IntVar(0, valmax, 'x6')
    x7 = solver.IntVar(0, valmax, 'x7')
    x8 = solver.IntVar(0, valmax, 'x8')

    print('Number of variables =', solver.NumVariables())

    solver.Add(2*x1+2*x5-3*x7>=0)
    solver.Add(2*x2+2*x6+6*x7>=0)
    solver.Add(4*x4-3*x7+4*x8>=0)
    solver.Add(x1+2*x2+2*x3-6*x8>=0)
    solver.Add(13*x3+3*x4-6*x5-4*x6>=0)
    solver.Add(3*x1+3*x2+6*x3+6*x4<=72)
    print('Number of constraints =', solver.NumConstraints())

    na=13
    nr=8
    nf=2
    solver.Maximize(2*na*x1+2*nr*x2+4*nf*x4+2*na*x5+2*nr*x6+(6*nr-3*na-3*nf)*x7+4*nf*x8)

    status = solver.Solve()

    if status == pywraplp.Solver.OPTIMAL:
        print('Solution:')
        print('Objective value =', solver.Objective().Value())
        print('x1 =', x1.solution_value())
        print('x2 =', x2.solution_value())
        print('x3 =', x3.solution_value())
        print('x4 =', x4.solution_value())
        print('x5 =', x5.solution_value())
        print('x6 =', x6.solution_value())
        print('x7 =', x7.solution_value())
        print('x8 =', x8.solution_value())
        xsol=[x1.solution_value(),x2.solution_value(),x3.solution_value(),x4.solution_value(),x5.solution_value(),x6.solution_value(),x7.solution_value(),x8.solution_value()]
        print('valeurs = ',statepath(xsol,[720.,0.,0.,0.,0.,0.]))
        print('valeurs = P,A,H,L,F,R')
    else:
        print('The problem does not have an optimal solution.')

    print('\nAdvanced usage:')
    print('Problem solved in %f milliseconds' % solver.wall_time())
    print('Problem solved in %d iterations' % solver.iterations())
    print('Problem solved in %d branch-and-bound nodes' % solver.nodes())
    return(xsol)


if __name__ == '__main__':
    main()

Number of variables = 8
Number of constraints = 6
Solution:
Objective value = 723.0
x1 = 0.0
x2 = 0.0
x3 = 12.0
x4 = 0.0
x5 = 26.0
x6 = 0.0
x7 = 5.0
x8 = 4.0
valeurs =  [  0. 370.   0.   0.  10. 300.]
valeurs = P,A,H,L,F,R

Advanced usage:
Problem solved in 7.000000 milliseconds
Problem solved in 0 iterations
Problem solved in 0 branch-and-bound nodes

-
Edité par Transistored 29 juillet 2020 à 10:35:01