ajout d'une image

[desynchronisation-controle.git] / exp_controle_asynchrone / simulMWSN.py~

diff --git a/exp_controle_asynchrone/simulMWSN.py~ b/exp_controle_asynchrone/simulMWSN.py~
index 171c682e28d84e56fdb608a54ebf20986e7acc32..863b2ae5bad14fea9dd48a0d9af268cacdb69da5 100644 (file)
--- a/exp_controle_asynchrone/simulMWSN.py~
+++ b/exp_controle_asynchrone/simulMWSN.py~
@@ -5,8 +5,11 @@ import numpy as np
 import pylab as pb
 from itertools import *
 from scipy import optimize as opt
 import pylab as pb
 from itertools import *
 from scipy import optimize as opt

+from copy import deepcopy 
+import sys as sy

 
 

-error  = 0.01
+
+error  = 0.1

 epsilon = 1E-10
 vrate = 0.8
 p = 0.7
 epsilon = 1E-10
 vrate = 0.8
 p = 0.7


@@ -16,7 +19,8 @@ nbiter = 1000
 POS = 1
 POS_NUL = 2
 POSINF1 = 3
 POS = 1
 POS_NUL = 2
 POSINF1 = 3

-
+init = []
+fichier_init="config_initiale_default.txt"

 
 
 
 
 
 


@@ -103,6 +107,7 @@ alpha = 0.5
 beta = 1.3E-8
 gamma = 55.54
 delta = 0.2
 beta = 1.3E-8
 gamma = 55.54
 delta = 0.2

+zeta = 0.1

 amplifieur = 1
 sigma2 = 3500
 Bi = 5
 amplifieur = 1
 sigma2 = 3500
 Bi = 5


@@ -215,18 +220,23 @@ def maj_theta(k):
 
 
 
 
 
 

-def maj(k,maj_theta):
+def maj(k,maj_theta,mxg,idxexp):

     # initialisation des operateurs lagrangiens
     # initialisation des operateurs lagrangiens

-    global u, v, la, w, theta , q,  Ps, Rh,  eta, x, valeurFonctionDuale 
-
-
+    global u, v, la, w, theta , q,  Ps, Rh,  eta, x, valeurFonctionDuale
+    
+    smax = 0

     #print "iteration",k
     
     up = {}        
     for h in V:
         for i in N:
     #print "iteration",k
     
     up = {}        
     for h in V:
         for i in N:

-            if not ASYNC or  random() < taux_perte:
-                up[(h,i)] = u[(h,i)]-theta*(eta[(h,i)]-sum([a[i][l]*x[(h,l)] for l in L]))
+            if not ASYNC or  random() < taux_succes:
+                s = eta[(h,i)]-sum([a[i][l]*x[(h,l)] for l in L])
+                if abs(s) > mxg :
+                    print "ds calcul u",abs(s),idxexp
+                    mxg = abs(s)                 
+                smax = max(smax,abs(s))
+                up[(h,i)] = u[(h,i)]-theta*s

             else :
                 up[(h,i)] = u[(h,i)]
 
             else :
                 up[(h,i)] = u[(h,i)]
 


@@ -234,10 +244,13 @@ def maj(k,maj_theta):
     
     vp = {}
     for h in V:
     
     vp = {}
     for h in V:

-        #print "vp",gamma*mt.pow(Ps[h],float(1)/3)
-        #vp[h] = max(0,v[h]-theta*(Rh[h]- mt.log(float(sigma2)/D)/(gamma*mt.pow(Ps[h],float(1)/3))))
-        if not ASYNC or  random() < taux_perte:
-            vp[h] = max(0,v[h]-theta*(Rh[h]- mt.log(float(sigma2)/D)/(gamma*mt.pow(Ps[h],float(1)/3))))
+        if not ASYNC or  random() < taux_succes:
+            s = Rh[h]- mt.log(float(sigma2)/D)/(gamma*mt.pow(Ps[h],float(1)/3))
+            if abs(s) > mxg :
+                print "ds calcul v",abs(s),idxexp
+                mxg = abs(s) 
+            smax = max(smax,abs(s))
+            vp[h] = max(0,v[h]-theta*s)

         else :
             vp[h] = v[h]
 
         else :
             vp[h] = v[h]
 


@@ -245,8 +258,13 @@ def maj(k,maj_theta):
 
     lap={}
     for i in N:
 
     lap={}
     for i in N:

-        if not ASYNC or  random() < taux_perte:
-            resla = la[i]-theta*((q[i]*Bi -sum([aplus[i][l]*cs[l]*sum([x[(h,l)] for h in V]) for l in L])-cr*sum([amoins[i][l]*sum([x[(h,l)] for h in V]) for l in L])-psi(Ps,i)))
+        if not ASYNC or  random() < taux_succes:
+            s = q[i]*Bi -sum([aplus[i][l]*cs[l]*sum([x[(h,l)] for h in V]) for l in L])-cr*sum([amoins[i][l]*sum([x[(h,l)] for h in V]) for l in L])-psi(Ps,i)
+            if abs(s) > mxg :
+                print "ds calcul la",abs(s),idxexp,i
+                mxg = abs(s) 
+            smax = max(smax,abs(s))
+            resla = la[i]-theta*s

             lap[i] = max(0,resla) 
         else :
             lap[i] = la[i]
             lap[i] = max(0,resla) 
         else :
             lap[i] = la[i]


@@ -256,8 +274,13 @@ def maj(k,maj_theta):
                         
     wp={}
     for l in L:
                         
     wp={}
     for l in L:

-        if not ASYNC or  random() < taux_perte:
-            wp[l] = w[l] + theta*(sum([a[i][l]*q[i] for i in N])) 
+        if not ASYNC or  random() < taux_succes:
+            s = sum([a[i][l]*q[i] for i in N])
+            if abs(s) > mxg :
+                print "ds calcul w",abs(s),idxexp
+                mxg = abs(s) 
+            smax = max(smax,abs(s))
+            wp[l] = w[l] + theta*s 

         else : 
             wp[l] = w[l]
 
         else : 
             wp[l] = w[l]
 


@@ -265,8 +288,6 @@ def maj(k,maj_theta):
 
     thetap = maj_theta(k)
 
 
     thetap = maj_theta(k)
 

-    #print "k",k,"maj theta", thetap
-

 
     qp={}
     def f_q(qi,i):
 
     qp={}
     def f_q(qi,i):


@@ -274,7 +295,7 @@ def maj(k,maj_theta):
         return qi*qi + qi*fa    
 
     for i in N:
         return qi*qi + qi*fa    
 
     for i in N:

-        if not ASYNC or  random() < taux_perte:
+        if not ASYNC or  random() < taux_succes:

             c = -float(sum([a[i][l]*w[l] for l in L]) - la[i]*Bi)/(2*amplifieur)
             rep = epsilon if c <= 0 else c
             qp[i] = rep
             c = -float(sum([a[i][l]*w[l] for l in L]) - la[i]*Bi)/(2*amplifieur)
             rep = epsilon if c <= 0 else c
             qp[i] = rep


@@ -291,10 +312,19 @@ def maj(k,maj_theta):
         #print "ds f_ps",psh, v[h]* mt.log(float(sigma2)/D)/(gamma*((psh**2)**(float(1)/3))) +la[h]*psh
         return v[h]* mt.log(float(sigma2)/D)/(gamma*mt.pow(float(2)/3)) +la[h]*psh
     for h in V:
         #print "ds f_ps",psh, v[h]* mt.log(float(sigma2)/D)/(gamma*((psh**2)**(float(1)/3))) +la[h]*psh
         return v[h]* mt.log(float(sigma2)/D)/(gamma*mt.pow(float(2)/3)) +la[h]*psh
     for h in V:

-        if not ASYNC or  random() < taux_perte:
+        if not ASYNC or  random() < taux_succes:
+            

             lah = 0.05 if la[h] == 0 else  la[h]
             lah = 0.05 if la[h] == 0 else  la[h]

-            rep = (float(2*v[h]*mt.log(float(sigma2)/D))/mt.pow(3*gamma*lah,float(3)/5))
+            rep = mt.pow(float(2*v[h]*mt.log(float(sigma2)/D))/(3*gamma*lah),float(3)/5)

             Psp[h] = epsilon if rep <= 0 else rep
             Psp[h] = epsilon if rep <= 0 else rep

+            """
+            t= float(-3*la[h]+mt.sqrt(9*la[h]**2+64*zeta*v[h]*mt.log(float(sigma2)/D)/gamma))/float(16*zeta)
+            #print t
+            rep = mt.pow(t,float(3)/5)
+            Psp[h]=rep
+            """
+############reprendre ici
+

         else :
             Psp[h] = Ps[h]
 
         else :
             Psp[h] = Ps[h]
 


@@ -313,7 +343,7 @@ def maj(k,maj_theta):
         return delta*rh*rh-v[h]*rh-sum([u[(h,i)]*eta[(h,i)] for i in N])
 
     for h in V:
         return delta*rh*rh-v[h]*rh-sum([u[(h,i)]*eta[(h,i)] for i in N])
 
     for h in V:

-        if not ASYNC or  random() < taux_perte:
+        if not ASYNC or  random() < taux_succes:

             rep = float(v[h])/(2*delta)
             Rhp[h] = 0 if rep < 0 else rep
         else : 
             rep = float(v[h])/(2*delta)
             Rhp[h] = 0 if rep < 0 else rep
         else : 


@@ -327,7 +357,7 @@ def maj(k,maj_theta):
     
     for h in V:
         for l in L:
     
     for h in V:
         for l in L:

-            if not ASYNC or  random() < taux_perte:
+            if not ASYNC or  random() < taux_succes:

                 rep = -float(cs[l]*sum([la[i]*aplus[i][l] for i in N]) +cr*sum([la[i]*amoins[i][l] for i in N])+sum([u[(h,i)]*a[i][l] for i in N]))/(2*delta)
                 xp[(h,l)] = 0 if rep < 0 else rep 
             else :
                 rep = -float(cs[l]*sum([la[i]*aplus[i][l] for i in N]) +cr*sum([la[i]*amoins[i][l] for i in N])+sum([u[(h,i)]*a[i][l] for i in N]))/(2*delta)
                 xp[(h,l)] = 0 if rep < 0 else rep 
             else :


@@ -351,7 +381,7 @@ def maj(k,maj_theta):
 
     arret = abs(valeurFonctionDuale-valeurFonctionDualep) < error
 
 
     arret = abs(valeurFonctionDuale-valeurFonctionDualep) < error
 

-    return (up,vp,lap,wp,thetap,etap,qp,Psp,Rhp,xp,valeurFonctionDualep,arret)
+    return (up,vp,lap,wp,thetap,etap,qp,Psp,Rhp,xp,valeurFonctionDualep,arret,mxg,smax)

 
 
 
 
 
 


@@ -373,48 +403,32 @@ print "x",x
 
 
 def initialisation():
 
 
 def initialisation():

-    global u, v, la, w, theta , q,  Ps, Rh,  eta, x, valeurFonctionDuale 
+    global u, v, la, w, theta , q,  Ps, Rh,  eta, x,init 

 
     theta = omega
 
 
     theta = omega
 

-
-    # TODO Ahmed : initialiser avec qi >0 et somme (ail.qi) = 0 Ei/i = Ti 
-    #q = [0 for i in N]

     q = {}
     for i in N :
         q[i] = 0.15 + random()*0.05
 
 
     q = {}
     for i in N :
         q[i] = 0.15 + random()*0.05
 
 

-
-
-
-    #TODO Ahmed doit donner les valeurs initiales pour Ps

     Ps= {}
     for h in V :
         Ps[h] =  0.2+random()*0.3  
 
     Ps= {}
     for h in V :
         Ps[h] =  0.2+random()*0.3  
 

-
-
-
-    #TODO Ahmed doit donner les valeurs initiales pour Rh

     Rh = {}
     for vi in V:
         Rh[vi] = 0.1 + random()*0.1
 
     Rh = {}
     for vi in V:
         Rh[vi] = 0.1 + random()*0.1
 

-
-

     eta = {}
     for h in V :
         for i in N:
             eta[(h,i)]= etahi(h,i,Rh)
 
     eta = {}
     for h in V :
         for i in N:
             eta[(h,i)]= etahi(h,i,Rh)
 

-
-
-    # Ok pour 0 

     x={}
     for h in V :
         for l in L:
     x={}
     for h in V :
         for l in L:

-            x[(h,l)]=0#random() 
+            x[(h,l)]=0

 
 
  
 
 
  


@@ -437,43 +451,94 @@ def initialisation():
     for l in L:
         w[l] =  random()
 
     for l in L:
         w[l] =  random()
 

+    init = [deepcopy(q),deepcopy(Ps),deepcopy(Rh),deepcopy(eta),
+            deepcopy(x),deepcopy(u),deepcopy(v),deepcopy(la),deepcopy(w)]

 
 
 
 

-def __evalue_maj_theta__():
+
+def initialisation_():
+    global u, v, la, w, theta , q,  Ps, Rh,  eta, x,init 
+    fd = open(fichier_init,"r")
+    l= fd.readline()
+    init_p = eval(l)
+    print init_p
+    theta = omega
+    (q,Ps,Rh,eta,x,u,v,la,w) = tuple([deepcopy(x) for x in init_p])
+    init = [deepcopy(q),deepcopy(Ps),deepcopy(Rh),deepcopy(eta),
+            deepcopy(x),deepcopy(u),deepcopy(v),deepcopy(la),deepcopy(w)]
+
+
+
+def __evalue_maj_theta__(nbexp,out=False):

     global u, v, la, w, theta , q,  Ps, Rh,  eta, x, valeurFonctionDuale 
     global u, v, la, w, theta , q,  Ps, Rh,  eta, x, valeurFonctionDuale 

-    nbexp = 50

     res = {}
     m = []
     res = {}
     m = []

+    itermax = 100000
+ 

     def __maj_theta(k):
     def __maj_theta(k):

+        mem = []

         om = omega/(mt.pow(k,0.75))
         return om
     for idxexp in range(nbexp):
         om = omega/(mt.pow(k,0.75))
         return om
     for idxexp in range(nbexp):

-        
-        initialisation()
+        mxg = 0
+        if not(out):
+            initialisation()
+        else :
+            initialisation_()
+            

         k = 1
         arret = False
         k = 1
         arret = False

-        while k < 1000000  and not arret : 
-            (u,v,la,w,theta,eta,q,Ps,Rh,x,valeurFonctionDuale,ar)=maj(k,__maj_theta)
-            errorq =  abs(min(q.values())-max(q.values()))
+        sm = 0
+        while k < itermax  and not arret : 
+            (u,v,la,w,theta,eta,q,Ps,Rh,x,valeurFonctionDuale,ar,mxg,smax)=maj(k,__maj_theta,mxg,idxexp)
+            errorq =  (max(q.values()) - min(q.values()))/min(q.values())

             arret = errorq <  error
             k+=1
             arret = errorq <  error
             k+=1

-            if k%5000 ==0 :
+            variation = "+" if smax > sm else "-"
+            sm = smax
+            print variation,
+            if k%100 ==0 :

                 print "k:",k,"erreur sur q", errorq, "et q:",q
                 print "k:",k,"erreur sur q", errorq, "et q:",q

+                print "maxg=", mxg
+                mem = [deepcopy(q),deepcopy(Ps),deepcopy(Rh),deepcopy(eta),
+                       deepcopy(x),deepcopy(u),deepcopy(v),deepcopy(la),deepcopy(w)]
+            if k%4500 == 0 :
+                print "#########\n",mem,"\#########\n"
+            if k%4600 == 0 :
+                print "#########\n",mem,"\#########\n"
+
+
+
+            if smax - sm  > 500:
+                print "variation trop grande"
+                print "init"
+                print init
+                sy.exit(0)
+        if k == itermax:
+            print "nbre d'iteration trop grand"
+            print "init"
+            print init
+            sy.exit(1)
+

         print "###############"
         print k
         print "###############"
         m += [k]
         print "###############"
         print k
         print "###############"
         m += [k]

+

     print (min(m),max(m),float(sum(m))/nbexp,m),m
 
 
 
 
 
     print (min(m),max(m),float(sum(m))/nbexp,m),m
 
 
 
 
 

-#ASYNC = False
-#__evalue_maj_theta()
-ASYNC = True
-taux_perte = 0.9
-__evalue_maj_theta__()
+ASYNC = False
+__evalue_maj_theta__(1,True)
+
+
+
+#ASYNC = True
+#taux_succes = 0.9
+#__evalue_maj_theta__()

 
 
 """
 
 
 """