Discussion :

[rcosentino]

Bonjour,

Je viens vers vous pour tenter d'optimiser un algorithme que j'ai écrit, qui fonctionne très bien, mais qui a le défaut d'être très lent.
Je suis absolument certain que l'on peut trouver une solution itérative plus efficace, mais comme cela sort de mon domaine de compétence, je m'en remet à vous.

Pour le fond de l'histoire, il s'agit de trouver la position d'équilibre d'une section en béton armé :


Pour faire simple, lorsqu'on soumet une poutre en béton armé à de la flexion (typiquement une poutre sur deux appuis supportant des charges descendantes), caractérisée par un moment de flexion noté M, on comprime la poutre en fibre supérieure et on la tend en fibre inférieure. La compression est équilibrée par la béton et la traction par les armatures qui sont noyées dans le béton.

Dans les cas complexes, la forme de la section peut être quelconque, la valeur de l'effort de compression dans le béton, que l'on notera Fc, n'est pas linéairement proportionnel à l'allongement de la fibre supérieure que l'on notera eps_sup.
De même, la valeur de l'effort de traction dans les armatures, que l'on notera Fs, n'est pas linéairement proportionnel à l'allongement de la fibre inférieure, que l'on notera eps_inf

On a :
f(eps_sup,eps_inf) = Fc
g(eps_sup,eps_inf) = Fs
h(eps_sup,eps_inf) = z

Les fonctions f, g et h sont non linéaires.

Lorsque l'équilibre est atteint, c'est à dire lorsqu'on a trouvé le bon couple de valeur eps_sup et eps_inf, les équations suivantes se vérifient :

Fc + Fs = 0 (avec Fs > 0 et Fs < 0)
M/z = Fc = -Fs (où M est une constante donnée, c'est le moment de flexion que l'on chercher à équilibrer)

Compte tenu des limites des matériaux, je connais les valeurs maximales et minimales admissibles pour eps_sup et eps_inf.
Je me donne un pas de calcul eps_step, suffisamment petit pour que le calcul converge systématiquement.

Mon algortihme, à ce jour, cherche la solution en recherchant le minimum de la fonction abs(Fc+Fs)/Fc + abs(M_calc-M)/M.

J'ai lu des choses intéressantes sur les algorithme de minimisation, mais je n'ai pas réussi à comprendre comme pratiquement programmer le calcul.

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
 
Constante : 
M
 
Variable :
delta = 10000000000000
M_calc = 0
eps_sup_candidat = 0
eps_inf_candidat = 0
 
 
Pour eps_sup = 0 à eps_sup_max par pas de eps_step
       Pour eps_inf = eps_min à 0 par pas de eps_step
                 Fc = f(eps_sup,eps_inf)
                 Fs = g(eps_sup,eps_inf)
                 z = h(eps_sup,eps_inf)
                 M_calc = Fc * z
 
                 Si (abs(Fc+Fs)/Fc + abs(M_calc-M)/M) < delta Alors
                         delta = abs(Fc+Fs)/Fc + abs(M_calc-M)/M
                         eps_sup_candidat = eps_sup
                         eps_inf_candidat = eps_inf
                 Fin Si
         Fin Pour
Fin Pour
 
Retrourner eps_sup_candidat et eps_inf_candidat

Par avance Merci !

[wiwaxia]

Bonjour,

... Je viens vers vous pour tenter d'optimiser un algorithme que j'ai écrit, qui fonctionne très bien, mais qui a le défaut d'être très lent.
Je suis absolument certain que l'on peut trouver une solution itérative plus efficace ...
Je me donne un pas de calcul eps_step, suffisamment petit pour que le calcul converge systématiquement ...

Tu es apparemment pris dans l'étau de deux contraintes opposées, qui te conduisent:
a) à choisir un pas (eps_step) aussi grand que possible, afin que le système parvienne rapidement à l'état d'équilibre recherché;
b) à éviter les valeurs trop élevées pour ce même pas, susceptibles de conduire à l'explosion numérique.

1°) Une solution consiste à envisager une augmentation progressive du pas à partir d'une valeur sûre (eps_step0) - celle-là même que tu emploies - au fur et à mesure de l'évolution du système; celle-ci doit cependant être paramétrée par un terme tendant vers zéro lorsque l'état d'équilibre est atteint.

De ce point de vue, le rapport sans dimension r = (Fc + Fs)/Fc pourrait très bien convenir compte tenu des informations données.

f(eps_sup,eps_inf) = Fc
g(eps_sup,eps_inf) = Fs
h(eps_sup,eps_inf) = z

Les fonctions f, g et h sont non linéaires.

Lorsque l'équilibre est atteint, c'est à dire lorsqu'on a trouvé le bon couple de valeur eps_sup et eps_inf, les équations suivantes se vérifient :

Fc + Fs = 0 (avec Fs > 0 et Fs < 0) (1)
M/z = Fc = -Fs (où M est une constante donnée, c'est le moment de flexion que l'on chercher à équilibrer) ...

(1): là, il y a eu lapsus: je suppose que l'on a d'après ce qui suit Fc > 0 .

2°) Il faut remarquer que (r) ne peut être rigoureusement nul - sauf cas numérique fortuit, parce que les calculs effectués se limitent à un nombre fini de chiffres (environ 16, le plus souvent); le rapport des forces Fs/Fc = r - 1 n'est connu qu'avec une incertitude au moins égale au "epsilon machine" du microprocesseur (Em, de l'ordre de 10^-16).

On pourra par exemple, après chaque centaine d'itérations effectuées, envisager une nouvelle valeur du pas à l'aide d'un algorithme du genre:

Code :

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

1
2
3
4
5
6
7
8
 
CONST eps_step1 = 1000 * eps_step0;
 
IF (r>1E-6) THEN eps_step:= eps_step0
            ELSE IF (r<1E-16) THEN eps_step:= eps_step1
                              ELSE BEGIN
                                     q:= Power(1E-6/r, 0.3); eps_step:= eps_step0 * q
                                   END;

qui permet d'augmenter continûment le pas d'un facteur 1000 lorsque le système s'approche de sa configuration d'équilibre.

[rcosentino]

Merci. J'avais pensé à ajuster le pas en fonction de la convergence. Mon problème, c'est que je code en VBA, et que la modification du pas n'est pas possible de façon dynamique dans une boucle for.

[tbc92]

Une boucle for peut toujours être remplacée par une boucle While :

Au lieu de :

Code :

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

1
2
3
for i = 0 to 1000 step 10
   //  traitement
next i

n
Tu peux faire :

Code :

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

1
2
3
4
5
6
7
i=0
step0 = 10
while i <= 1000
  //  traitement
  if ... then step0 = step0+1
  i = i+step0
end while

[rcosentino]

Merci.

J’ai divisé mon temps de calcul par 20. C’est parfait.

[Invité]

hello,

J'ai lu des choses intéressantes sur les algorithme de minimisation, mais je n'ai pas réussi à comprendre comme pratiquement programmer le calcul.

effectivement tu peux tenter de minimiser g(x) = abs((1))+abs((2))
(avec x le vecteur (variable) [epsinf, epssup]
où (1) = Fc-Fs et (2) = M/z-Fc

le truc c'est que prendre la valeur absolue c'est pas continu... donc tu peux prendre une fonction sym en 0 continue par ex
x^2
idem g(x) = (Fc-Fs)^2 + (M/z-Fc)^2
(et g reste minimale qd tes deux termes sous leur carré valent 0)

Ensuite tu peux appliquer une bête descente de gradient sur g (grosso faut calculer g'(x))
tu peux poser plutot

g(x) = 1/2 [ (Fc-Fs)^2 + (M/z-Fc)^2]

puisqu'en 1/2 changera rien à la convergence mais simplifie pour la dérivation

g'(x) = (Fc'(x)-Fs'(x)).*(Fc(x)-Fs(x)) + (-z'(x)/z^2(x)*M - Fc'(x)).*(M/z(x)-Fc(x))

a noter ici je fais un ecart de notation, la dériviation est partielle sur les variables respectives epsinf et epssup


SI tu sais calculer analytiquement Fs'(x) et z'(x) c'est tranquillou bilou tu regardes wiki gradient descent (et t'adaptes pr la forme vectorielle ou lieu de scalaire)
idem

SI l'expression de Fc et z est trop chaude tu peux qd même () t'en sortir en les calculant discrètement mais ya des truc un peu complexes donc le plus "serein" c'est tu calcules
g'(x) = (g(x+h) - g(x))./h avec h que tu prends bien petit (et de la forme 1e-8(1,1) par ex)

ds le code adapté ci-dessous de https://en.wikipedia.org/wiki/Gradient_descent, df(x,y) correspond à g'([epsinf, epssup])

Code :

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

1
2
3
4
5
6
7
8
9
10
11
 
    df(x,y):
        h=1e-8
        return [ (f(x+h)-f(x))/h , (f(y+h)-f(y))/h ]
 
    while (previous_step_size > precision) & (iters < max_iters):
        prev_x = cur_x
        prev_y = cur_y
        cur_x, cur_y -= gamma * df(prev_x, prev_y)
        previous_step_size = sqrt( (cur_x - prev_x)^2 + (cur_y - prev_y)^2)
        iters+=1

A s'assurer que tu converges vers là où tu veux (à priori initialiser (cur_x,cur_y)=(epsmin/2, epsmax/2) à l'air pas trop bête: tu commences déjà dans ton intervalle "solution" et tu n'as plus qu'à vérifier à la convergence si tu y es bien resté. Sinon tu n'avais pas d'optimum local dans ton intervalle solution et donc pas de solution respectant tes contraintes pour epsinf et epssup)