Problème Boucle dans Boucle

**Sniffle** · 10/06/2009, 14h53

Bonjour à vous,

J'ai besoin de votre aide pour une imbrication de boucles.

J'ai donc deux boucles imbriquées l'une dans l'autre de la façon suivante :

Boucle Géométrie
Boucle Stabilisation

Fin de boucle Stabilisation
Fin de boucle Géométrie

Pour vous expliquer le fonctionnement du programme : J'ai deux câbles qui ont une influence l'un sur l'autre. La position du câble du bas doit respecter une position objectif tandis que celle du haut dépend des efforts ponctuels nécessaires à l'atteinte de la position objectif du câble du bas.

Les câbles sont découpés en éléments et en noeuds, et il y a des forces ponctuelles tous les n noeuds.

Mon problème est dans la boucle géométrie. Car en effet, la boucle Stabilisation fonctionne quand je n'utilise qu'elle. Et quand je lui rajoute la boucle géométrie (Qui est une boucle qui recherche la solution pas à pas par ajout/retrait de Delta force à la valeur des forces ponctuelles pour obtenir une convergence et qui fonctionne dans le cas simple de recherche de solution à des équations basiques) le programme se bloque.

La blocage intervient de la manière suivante :

Le calcul de la première stabilisation se fait entièrement (30 calculs). Le programme passe donc à la vérification (Je suppose) de la réalisation ou non de l'objectif géométrie et ensuite, il effectue le premier calcul de la boucle stabilisation à plusieurs reprises.

Je vous donne en dessous les boucles en questions :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
def boucle_Newton_Raphson(etat,poi,pen):
 
    geometrie = False
    forme_objectif = 1.4
    compteur = 0
 
    while not(geometrie):
 
        compteur+=1
        fini = False
        compt = 0
 
        while not(fini):
            compt += 1
 
            ### calcul de K et des efforts internes
            (KK_glob,ff_glob) = calcule_K_f(etat,liste_elems)
 
            ### calcul des efforts externes
            # cas simple 100 N par noeud vers le bas
            f_ext = matrix(zeros((2*nb_noeuds,1)))
 
            for i in range(nb_noeuds):
                f_ext[2*i+1,0] = - carac_ml*dx
 
            ### Pendules
            force_actuelle=[0]             
            f_pendule = matrix(zeros((2*nb_noeuds,1)))    
            for i in range (nb_pendules):
                f_pendule[i*2*npen+1,0]= poi * carac_ml*ltot / nb_pendules - pen       
                force_actuelle.append(f_pendule[i*2*npen+1,0])
 
            #f_pendule = matrix(zeros((2*nb_noeuds,1)))
            #for i in range(nb_pendules):
                #f_pendule[i*2*npen+1,0] = - f_ext[2*i+1,0] * npen * poi - pen
 
            # for n in range(nb_pendules):
             #   f_pendule[2*((1+n)*int(nb_noeuds/(nb_pendules+1)))+1,0] = - f_ext[2*((1+n)*int(nb_noeuds/(nb_pendules+1)))+1,0] * nb_noeuds/(nb_pendules+1) * poi/2 - pen
 
 
            ### reduction du systeme
            (K_red,f_red) = systeme_reduit(KK_glob,ff_glob+f_ext+f_pendule)#+f_inverse)#f_surcharge)
 
            ### resolution
            u_compl = deplacement_complet(K_red,f_red)
 
            ### affichage
            print("Calcul "+str(compt)+" -- min(etat)="+str(min(etat)[0,0]))
 
            ### test de convergence
            if max(abs(u_compl))<1e-10:
                fini = True
 
            # limitation des deplacements si necessaire
            tol = 0.5
            if max(abs(u_compl)) > tol:
                u_compl *= tol/max(abs(u_compl))[0,0]
 
            #trace graphe
            #trace_graphe('b')
            #draw()
 
            ### mise a jour etat
            etat = etat + u_compl
            for i in range(nb_pendules):
                print(etat[i*2*npen+1,0])
 
    ##premiere boucle            
    for i in range(nb_pendules):
        if etat[i*2*npen+1,0] > (forme_objectif-0.1) and etat(i*2*npen+1,0) < (forme_objectif+0.1):
            geometrie = True
 
        else :
            geometrie = False
            fini = False
            for i in range (nb_pendules):
                print (i)
                if etat[i*2*npen+1,0] > forme_objectif:  
                    if etat[i*2*npen+1,0]-forme_objectif > 1:
                        force_actuelle[i]-=0.1
                    else:
                        if etat[i*2*npen+1,0]-forme_objectif > .3 :
                            force_actuelle[i]-=0.05
                        else:
                            force_actuelle[i]-=0.0001                
                else:
                    if etat[i*2*npen+1,0]-forme_objectif < 1:
                        force_actuelle[i]+=0.1
                    else:
                        if etat[i*2*npen+1,0]-forme_objectif < .3 :
                            force_actuelle[i]+=0.05
                        else:
                            force_actuelle[i]+=0.0001 
 
                print("1",compt, force_actuelle[i], etat[i*2*npen+1,0],0) 
 
    return etat

De plus, je voudrais que pour un câble, le calcul tienne compte des deux boucles car il y a une géométrie objectif pour lui et que pour le second câble on ne tienne compte que des forces ponctuelles obtenues pour le premier câble et qu'il n'y ait donc qu'un passage dans la boucle stabilisation et pas la boucle géométrie car il n'y a pas d'objectif de ce côté là.

Merci pour vos lumières,

Sniffle

**eyquem** · 17/06/2009, 13h57

Salut.

Grâce à ton effort louable de présenter le problème clairement, j’ai bien compris jusqu’à «Le blocage intervient de la manière suivante».
Ensuite je ne comprends plus:

- qu’est ce que tu appelles première stabilisation ? quelles sont les lignes du code concernées ?

- où est la vérification ?
est-ce
if etat[i*2*npen+1,0] > (forme_objectif-0.1) and etat(i*2*npen+1,0) < (forme_objectif+0.1): ?

- «premier calcul de la boucle stabilisation» : quelles sont les lignes de la boucle stabilisation ? ?

- pourquoi la dernière partie est-elle appelée «première boucle» ?

- la dernière phrase «De plus, je voudrais que pour un câble, le calcul tienne compte des deux boucles car il y a une géométrie objectif pour lui et que pour .....» m’est totalement incompréhensible

Toutefois:

1/ j’ai comme l’impression que la dernière partie, celle qui commence à «## premiere boucle» est mal indentée:

si la variable fini ne prend sa valeur False que dans cette partie, cette dernière doit être à l’intérieur de la boucle while not (fini) puisque la variable fini est le critère d’arrêt de cette boucle.

Il faudrait donc l’indenter 2 fois vers la droite de manière qu'il y ait cet alignement

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
 
while not(fini):
    compt += 1
    etc
    .......
    ##premiere boucle 
    for i in range(nb_pendules):

Mais je n’en suis pas très sûr parce que je ne comprends ce code que superficiellement.

2/ J’ai l’impression qu’il y a aussi un problème de conception algorithmique.

Je me demande si le code suivant ne correspond pas mieux à ce que tu cherches à faire. Là encore, je ne suis pas sûr du tout.
(la modif est en bas au niveau de geometrie = True)

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
def boucle_Newton_Raphson(etat,poi,pen):
 
    geometrie = False
    forme_objectif = 1.4
    compteur = 0
 
    while not(geometrie):
 
        compteur+=1
        fini = False
        compt = 0
 
        while not(fini):
            compt += 1
 
            ### calcul de K et des efforts internes
            (KK_glob,ff_glob) = calcule_K_f(etat,liste_elems)
 
            ### calcul des efforts externes
            # cas simple 100 N par noeud vers le bas
            f_ext = matrix(zeros((2*nb_noeuds,1)))
 
            for i in range(nb_noeuds):
                f_ext[2*i+1,0] = - carac_ml*dx
 
            ### Pendules
            force_actuelle=[0]             
            f_pendule = matrix(zeros((2*nb_noeuds,1)))    
            for i in range (nb_pendules):
                f_pendule[i*2*npen+1,0]= poi * carac_ml*ltot / nb_pendules - pen       
                force_actuelle.append(f_pendule[i*2*npen+1,0])
 
            #f_pendule = matrix(zeros((2*nb_noeuds,1)))
            #for i in range(nb_pendules):
                #f_pendule[i*2*npen+1,0] = - f_ext[2*i+1,0] * npen * poi - pen
 
            # for n in range(nb_pendules):
             #   f_pendule[2*((1+n)*int(nb_noeuds/(nb_pendules+1)))+1,0] = - f_ext[2*((1+n)*int(nb_noeuds/(nb_pendules+1)))+1,0] * nb_noeuds/(nb_pendules+1) * poi/2 - pen
 
 
            ### reduction du systeme
            (K_red,f_red) = systeme_reduit(KK_glob,ff_glob+f_ext+f_pendule)#+f_inverse)#f_surcharge)
 
            ### resolution
            u_compl = deplacement_complet(K_red,f_red)
 
            ### affichage
            print("Calcul "+str(compt)+" -- min(etat)="+str(min(etat)[0,0]))
 
            ### test de convergence
            if max(abs(u_compl))<1e-10:
                fini = True
 
            # limitation des deplacements si necessaire
            tol = 0.5
            if max(abs(u_compl)) > tol:
                u_compl *= tol/max(abs(u_compl))[0,0]
 
            #trace graphe
            #trace_graphe('b')
            #draw()
 
            ### mise a jour etat
            etat = etat + u_compl
            for i in range(nb_pendules):
                print(etat[i*2*npen+1,0])
 
            krit = [ False for i in xrange(nb_pendules) ]        
            for i in range(nb_pendules):
                if (forme_objectif-0.1) < etat(i*2*npen+1,0) < (forme_objectif+0.1):
                    krit[i] = True
            if all(krit):
                geometrie = True
 
            else :
                geometrie = False
                fini = False
                for i in range (nb_pendules):
                    print (i)
                    if etat[i*2*npen+1,0] > forme_objectif:  
                        if etat[i*2*npen+1,0]-forme_objectif > 1:
                            force_actuelle[i]-=0.1
                        else:
                            if etat[i*2*npen+1,0]-forme_objectif > .3 :
                                force_actuelle[i]-=0.05
                            else:
                                force_actuelle[i]-=0.0001                
                    else:
                        if etat[i*2*npen+1,0]-forme_objectif < 1:
                            force_actuelle[i]+=0.1
                        else:
                            if etat[i*2*npen+1,0]-forme_objectif < .3 :
                                force_actuelle[i]+=0.05
                            else:
                                force_actuelle[i]+=0.0001 
 
                    print("1",compt, force_actuelle[i], etat[i*2*npen+1,0],0) 
 
    return etat

Si c’est le cas, je modifierais bien ce code de la manière suivante:
j'elimine la variable geometrie, je mets while 1, je simplifie encore la ligne de test, même si elle te paraît plus compliquée.

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
def boucle_Newton_Raphson(etat,poi,pen):
 
    forme_objectif = 1.4
    compteur = 0
 
    while 1:
 
        compteur+=1
        fini = False
        compt = 0
 
        while not(fini):
            compt += 1
 
            ### calcul de K et des efforts internes
            (KK_glob,ff_glob) = calcule_K_f(etat,liste_elems)
 
            ### calcul des efforts externes
            # cas simple 100 N par noeud vers le bas
            f_ext = matrix(zeros((2*nb_noeuds,1)))
 
            for i in range(nb_noeuds):
                f_ext[2*i+1,0] = - carac_ml*dx
 
            ### Pendules
            force_actuelle=[0]             
            f_pendule = matrix(zeros((2*nb_noeuds,1)))    
            for i in range (nb_pendules):
                f_pendule[i*2*npen+1,0]= poi * carac_ml*ltot / nb_pendules - pen       
                force_actuelle.append(f_pendule[i*2*npen+1,0])
 
            #f_pendule = matrix(zeros((2*nb_noeuds,1)))
            #for i in range(nb_pendules):
                #f_pendule[i*2*npen+1,0] = - f_ext[2*i+1,0] * npen * poi - pen
 
            # for n in range(nb_pendules):
             #   f_pendule[2*((1+n)*int(nb_noeuds/(nb_pendules+1)))+1,0] = - f_ext[2*((1+n)*int(nb_noeuds/(nb_pendules+1)))+1,0] * nb_noeuds/(nb_pendules+1) * poi/2 - pen
 
 
            ### reduction du systeme
            (K_red,f_red) = systeme_reduit(KK_glob,ff_glob+f_ext+f_pendule)#+f_inverse)#f_surcharge)
 
            ### resolution
            u_compl = deplacement_complet(K_red,f_red)
 
            ### affichage
            print("Calcul "+str(compt)+" -- min(etat)="+str(min(etat)[0,0]))
 
            ### test de convergence
            if max(abs(u_compl))<1e-10:
                fini = True
 
            # limitation des deplacements si necessaire
            tol = 0.5
            if max(abs(u_compl)) > tol:
                u_compl *= tol/max(abs(u_compl))[0,0]
 
            #trace graphe
            #trace_graphe('b')
            #draw()
 
            ### mise a jour etat
            etat = etat + u_compl
            for i in range(nb_pendules):
                print(etat[i*2*npen+1,0])
 
 
            if all([(forme_objectif-0.1 < etat(i*2*npen+1,0) < forme_objectif+0.1) for for i in range(nb_pendules)]:
                return etat
 
            else :
                fini = False
                for i in range (nb_pendules):
                    print (i)
                    if etat[i*2*npen+1,0] > forme_objectif:  
                        if etat[i*2*npen+1,0]-forme_objectif > 1:
                            force_actuelle[i]-=0.1
                        else:
                            if etat[i*2*npen+1,0]-forme_objectif > .3 :
                                force_actuelle[i]-=0.05
                            else:
                                force_actuelle[i]-=0.0001                
                    else:
                        if etat[i*2*npen+1,0]-forme_objectif < 1:
                            force_actuelle[i]+=0.1
                        else:
                            if etat[i*2*npen+1,0]-forme_objectif < .3 :
                                force_actuelle[i]+=0.05
                            else:
                                force_actuelle[i]+=0.0001 
 
                    print("1",compt, force_actuelle[i], etat[i*2*npen+1,0],0)

Si le critère qui execute geometrie = True dans ton code et donne all([etc]) ==True dans le mien est vérifié, l’instruction return fait sortir immediatement de toute la fonction. Il n’y a pas besoin de break, ou autre conditionnalité sur une variable.

3/ Il se pourrait que le blocage soit dû au fait que la dernière partie n’est pas correctement écrite. Là j’en suis sûr: il y a une erreur d’écriture. Mais ça ne veut pas dire que c’est le seul problème dans ton code

Si etat[i*2*npen+1,0] > forme_objectif n’est pas vérifié, alors c’est que
etat[i*2*npen+1,0] < forme_objectif

Donc si on teste ensuite
etat[i*2*npen+1,0]-forme_objectif < 1
c’est à dire si on teste
etat[i*2*npen+1,0]< forme_objectif + 1

c’est sûr que ce test donnera toujours True.

Par conséquent, quand etat[i*2*npen+1,0] < forme_objectif, la correction de force_actuelle[i] est toujours faite par +=0.1 : ça ne peut jamais s’approcher correctement de la solution

Ci dessous le code corrigé à ce niveau et réécrit plus clairement avec elif:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
          else :
                geometrie = False
                fini = False
                for i in range (nb_pendules):
                    print (i)
 
                    if etat[i*2*npen+1,0] > forme_objectif:  
                        if etat[i*2*npen+1,0]-forme_objectif > 1:
                            force_actuelle[i]-=0.1
                        elif etat[i*2*npen+1,0]-forme_objectif > .3 :
                            force_actuelle[i]-=0.05
                        else:
                            force_actuelle[i]-=0.0001
 
                    else: # etat[i*2*npen+1,0] < forme_objectif
                        if forme_objectif - etat[i*2*npen+1,0] > 1:
                            force_actuelle[i]+=0.1
                        elif etat[i*2*npen+1,0]forme_objectif - etat[i*2*npen+1,0] > .3 :
                            force_actuelle[i]+=0.05
                        else:
                            force_actuelle[i]+=0.0001 
 
                    print("1",compt, force_actuelle[i], etat[i*2*npen+1,0],0)

**Sniffle** · 23/06/2009, 16h32

Merci pour la réponse, je n'en attentais plus.

Dans la première version de mon programme, la seule boucle qui était présente était une boucle de recherche de la position d'équilibre qui avait comme condition d'arrêt : u_compl < 1e-10. C'est elle que j'appelle la première stabilisation.

c_compl est le déplacement du câble.

La condition de convergence est donc une condition sur le déplacement.

Ensuite, pour améliorer le programme, j'ai voulu ajouter une seconde boucle qui permettait de réitérer le calcul tant que la condition de convergence que tu as mise en évidence (if etat[i*2*npen+1,0] > (forme_objectif-0.1) and etat(i*2*npen+1,0) < (forme_objectif+0.1) n'était pas réalisée.

Ceci est une condition sur la géométrie du câble du bas.

L'appellation «premiere boucle» est le vestige de plusieurs copiés/collés.

Sinon, pour la phrase que tu ne parviens pas à comprendre, je vais essayer d'être plus clair.

J'ai donc deux câbles :

Un câble porteur qui n'a aucune contrainte de géométrie.

Un câble de contact qui doit avoir une forme précise.

Par conséquent, le calcul pour le câble de contact devra se faire en tenant compte des deux boucles While (stabilisation «fini», géométrie)

Alors que le calcul pour le câble porteur, n'a besoin que de la boucle While de stabilisation «Fini».

Mais en y réfléchissant, la réponse est toute bête pour ce problème. Il suffit de mettre des conditions de convergence pour la boucle géométrie, très faciles à réaliser.

Sinon, je vais analyser ce que tu m'as dit. J'aurais du jeter un oeil plus tôt car la fin de mon projet approche. J'ai trouvé une autre solution qui fonctionne mais celle là me permettrait d'aller plus loin donc encore merci de ta réponse.

Problème Boucle dans Boucle

Calcul scientifique Python

Discussions similaires

Partager

Partager