exp_controle_asynchrone/simulMWSNASYNC.py

   1 import math as mt
   2 from random import *
   3 import networkx as nx
   4 import numpy as np
   5 import pylab as pb
   6 from itertools import *
   7 from scipy import optimize as opt
   8
   9 error  = 1E-5
  10 epsilon = 1E-10
  11 vrate = 0.8
  12 p = 0.7
  13 coteCarre = 50
  14 distanceEmmissionMax = 30
  15 nbiter = 1000
  16 POS = 1
  17 POS_NUL = 2
  18 POSINF1 = 3
  19
  20 ASYNC = False
  21 taux_perte = 0.99
  22
  23 # construction du  graphe
  24 n = 10
  25 lg = [(0, 1, 22.004323820359151), (0, 2, 28.750632705280324), (0, 3, 29.68069293796183), (0, 4, 8.547146256271331), (0, 5, 28.079672647730469), (0, 7, 23.017867703525138), (0, 8, 6.1268526078857208), (0, 9, 24.573433868296771), (1, 0, 22.004323820359151), (1, 2, 18.807277287689722), (1, 3, 18.982897767602783), (1, 4, 16.848855991756174), (1, 5, 17.042671653231526), (1, 6, 16.410544777532913), (1, 7, 25.598808236367063), (1, 8, 20.175759189503321), (1, 9, 12.843763853990932), (2, 0, 28.750632705280324), (2, 1, 18.807277287689722), (2, 3, 1.0957062702237066), (2, 4, 29.159997765424084), (2, 5, 1.8557839901886808), (2, 6, 23.122260476726876), (2, 9, 6.052562826627808), (3, 0, 29.68069293796183), (3, 1, 18.982897767602783), (3, 2, 1.0957062702237066), (3, 4, 29.884008054261855), (3, 5, 1.9922790489539697), (3, 6, 22.479228556182363), (3, 9, 6.4359869969688059), (4, 0, 8.547146256271331), (4, 1, 16.848855991756174), (4, 2, 29.159997765424084), (4, 3, 29.884008054261855), (4, 5, 28.006189408396626), (4, 7, 15.774839848636024), (4, 8, 3.6206480052249144), (4, 9, 23.804744370383144), (5, 0, 28.079672647730469), (5, 1, 17.042671653231526), (5, 2, 1.8557839901886808), (5, 3, 1.9922790489539697), (5, 4, 28.006189408396626), (5, 6, 21.492976178079076), (5, 8, 29.977996181215822), (5, 9, 4.4452006140146185), (6, 1, 16.410544777532913), (6, 2, 23.122260476726876), (6, 3, 22.479228556182363), (6, 5, 21.492976178079076), (6, 9, 20.04488615603487), (7, 0, 23.017867703525138), (7, 1, 25.598808236367063), (7, 4, 15.774839848636024), (7, 8, 16.915923579829141), (8, 0, 6.1268526078857208), (8, 1, 20.175759189503321), (8, 4, 3.6206480052249144), (8, 5, 29.977996181215822), (8, 7, 16.915923579829141), (8, 9, 25.962918470558208), (9, 0, 24.573433868296771), (9, 1, 12.843763853990932), (9, 2, 6.052562826627808), (9, 3, 6.4359869969688059), (9, 4, 23.804744370383144), (9, 5, 4.4452006140146185), (9, 6, 20.04488615603487), (9, 8, 25.962918470558208)]
  26
  27 #n=3
  28 #lg= [(0,1,23),(1,0,15),(1,2,45)]
  29 sink = n-1
  30
  31
  32 def genereGraph():
  33     test = False
  34     while not test :
  35         G = nx.DiGraph()
  36         l = [(random()*coteCarre, random()*coteCarre) for _  in range(n)]
  37         for io in range(len(l)) :
  38             for ie in range(len(l)) :
  39                 if ie != io  :
  40                     dist = mt.sqrt((l[io][0]-l[ie][0])**2 + (l[io][1]-l[ie][1])**2)
  41                     if dist <= distanceEmmissionMax :
  42                         G.add_edge(io,ie,weight=dist)
  43                         G.add_edge(ie,io,weight=dist)
  44         test = not(any([ not(nx.has_path(G,o,sink)) for o in  G.nodes() if sink in G.nodes() and o != sink]))
  45     return G
  46
  47 G = nx.DiGraph()
  48 G.add_weighted_edges_from(lg)
  49 #print nx.is_strongly_connected(G)
  50
  51
  52 #nx.draw(G)
  53 #pb.show()
  54
  55
  56
  57
  58
  59
  60 #print G.edges(data=True)
  61 #TODO afficher le graphe et etre sur qu'il est connexe
  62
  63 N = G.nodes()
  64
  65
  66 #V = list(set(sample(N,int(len(N)*vrate)))-set([sink]))
  67 V = list(set(N)-set([sink]))
  68
  69
  70 source = V
  71 #print "source",source
  72
  73
  74 L = range(len(G.edges()))
  75 d = [di['weight'] for (_,_,di) in G.edges(data=True)]
  76
  77
  78
  79 def ail(i,l):
  80     assert  l in L, " pb de dimennsion de l: "+str(l)+ " "+ str(L)
  81     r = 0
  82     (o,e) = G.edges()[l]
  83     if o == i :
  84         r = 1
  85     elif e == i :
  86         r = -1
  87     return r
  88
  89
  90
  91 # Constantes
  92 a = [[ail(i,l) for l in L ] for i in xrange(n)]
  93 aplus  = [[1 if  ail(i,l) > 0 else 0  for l in L ] for i in xrange(n)]
  94 amoins  = [[1 if  ail(i,l) < 0 else 0  for l in L ] for i in xrange(n)]
  95
  96
  97
  98
  99
 100 alpha = 0.5
 101 beta = 1.3E-8
 102 gamma = 55.54
 103 delta = 0.2
 104 amplifieur = 1
 105 sigma2 = 3500
 106 Bi = 5
 107 omega = 0.15
 108 D = 100
 109 path_loss_exp = 4
 110 cr = 0.5
 111 cs = [alpha + beta*(di**path_loss_exp) for di in d]
 112 #print "cs",cs
 113
 114
 115 #TODO
 116 def etahi(h,i,R):
 117     ret = 0
 118     if i == h :
 119         ret = R[h]
 120     elif i == sink  :
 121         ret = -R[h]
 122     return ret
 123
 124
 125 def psi(Ps,i):
 126     if i not in Ps:
 127         return 0
 128     else :
 129         return Ps[i]
 130
 131 # separer le calcul des operateurs lagrangien pour avoir un pas d'iteration de decallage
 132
 133
 134 def cmpPair(x1,x2):
 135     return cmp(x1[1],x2[1])
 136
 137 def aminp(l):
 138     l.sort(cmp=cmpPair)
 139     return l[0][0]
 140
 141
 142
 143 def distance(d1,d2):
 144     return mt.sqrt(sum([(d1[t]-d2[t])**2 for t in d1]))
 145
 146
 147 def AfficheVariation (up,vp,lap,wp,thetap,etap,qp,Psp,Rhp,xp,valeurFonctionDualep):
 148     global u, v, la, w, theta , q, Ps, Rh, eta, x,valeurFonctionDuale
 149
 150     print "du=",distance(u,up),
 151     print "dv=",distance(v,vp),
 152     print "dlambda=",distance(la,lap),
 153     print "dw=",distance(w,wp),
 154     print "dtheta=",abs(theta-thetap),
 155     print "deta=",distance(eta,etap),
 156     print "dq=",distance(q,qp),
 157     print "dPs=",distance(Ps,Psp),
 158     print "dRh=",distance(Rh,Rhp),
 159     print "dx=",distance(x,xp),
 160     print "dL=",abs(valeurFonctionDuale-valeurFonctionDualep),"\n"
 161
 162
 163 valeurFonctionDuale = 0
 164
 165 def entre0et1(x):
 166     r = x if (x >0 and x <= 1) else -1
 167     #print "ds entre0et1 x et r", x,r
 168     return r
 169
 170
 171 def xpos(x):
 172     r = x if x >0 else -1
 173     #print "ds xpos x et r", x,r
 174     return r
 175
 176 def xposounul(x):
 177     return x if x >= 0 else -1
 178
 179 #f_q,q[i],POS,[i]
 180 def armin(f,xini,xr,args):
 181     #xpos = lambda x : x if x > 0 else -1
 182     r = 0
 183     if xr == POS :
 184         #print "strictement pos"
 185         #print "parametres passes a cobyla",xini,xpos,args,"consargs=(),rhoend=1E-5,iprint=0,maxfun=1000"
 186         r= opt.fmin_cobyla(f,xini,cons=[xpos],args=args,consargs=(),rhoend=1E-3,iprint=0,maxfun=nbiter)
 187         #print "le min str pos est",r
 188         #print "le min pos   est", r,"avec xini",xini
 189     elif xr == POS_NUL :
 190         #print "pos ou nul"
 191         r = opt.fmin_cobyla(f,xini,[xposounul],args,consargs=(),rhoend=1E-3,iprint=0,maxfun=nbiter)
 192         # print "le min pos est", r
 193         #print "le min pos null  est", r,"avec xini",xini
 194     elif xr == POSINF1:
 195         r = opt.fmin_cobyla(f,xini,[entre0et1],args,consargs=(),rhoend=1E-3,iprint=0,maxfun=nbiter)
 196         #print "le min pos inf 1 est", r,"avec xini",xini
 197
 198
 199     return r
 200
 201
 202
 203
 204
 205
 206 def maj_theta(k):
 207     return omega/(mt.pow(k,0.5))
 208
 209
 210
 211 def maj_theta(k):
 212     return omega/mt.sqrt(k)
 213
 214
 215
 216
 217
 218 def maj(k,maj_theta):
 219     # initialisation des operateurs lagrangiens
 220     global u, v, la, w, theta , q,  Ps, Rh,  eta, x, valeurFonctionDuale
 221
 222
 223     #print "iteration",k
 224
 225     up = {}
 226     for h in V:
 227         for i in N:
 228             if not ASYNC or  random() < taux_perte:
 229                 up[(h,i)] = u[(h,i)]-theta*(eta[(h,i)]-sum([a[i][l]*x[(h,l)] for l in L]))
 230             else :
 231                 up[(h,i)] = u[(h,i)]
 232
 233
 234
 235
 236     vp = {}
 237     for h in V:
 238         #print "vp",gamma*mt.pow(Ps[h],float(1)/3)
 239         if not ASYNC or  random() < taux_perte:
 240             vp[h] = max(0,v[h]-theta*(Rh[h]- mt.log(float(sigma2)/D)/(gamma*mt.pow(Ps[h],float(1)/3))))
 241         else :
 242             vp[h] = v[h]
 243
 244
 245     lap={}
 246     for i in N:
 247         if not ASYNC or  random() < taux_perte:
 248             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)))
 249             lap[i] = max(0,resla)
 250         else :
 251             lap[i] = la[i]
 252
 253     wp={}
 254     for l in L:
 255         if not ASYNC or  random() < taux_perte:
 256             wp[l] = w[l] + theta*(sum([a[i][l]*q[i] for i in N]))
 257         else :
 258             wp[l] = w[l]
 259
 260     thetap = maj_theta(k)
 261
 262
 263     qp={}
 264     def f_q(qi,i):
 265         fa = sum([a[i][l]*w[l] for l in L]) - la[i]*Bi
 266         return qi*qi + qi*fa
 267
 268     for i in N:
 269         if not ASYNC or  random() < taux_perte:
 270             c = -float(sum([a[i][l]*w[l] for l in L]) - la[i]*Bi)/(2*amplifieur)
 271             rep = epsilon if c <= 0 else c
 272             qp[i] = rep
 273         else :
 274             qp[i] = q[i]
 275
 276
 277
 278
 279
 280
 281
 282     Psp={}
 283     def f_Ps(psh,h):
 284         return v[h]* mt.log(float(sigma2)/D)/(gamma*mt.pow(float(2)/3)) +la[h]*psh
 285     for h in V:
 286         if not ASYNC or  random() < taux_perte:
 287             lah = 0.05 if la[h] == 0 else  la[h]
 288             rep = (float(2*v[h]*mt.log(float(sigma2)/D))/mt.pow(3*gamma*lah,float(3)/5))
 289             Psp[h] = epsilon if rep <= 0 else rep
 290         else :
 291             Psp[h] = Ps[h]
 292
 293
 294
 295
 296     etap={}
 297     for h in V:
 298         for i in N :
 299             etap[(h,i)] = etahi(h,i,Rh)
 300
 301
 302
 303     Rhp={}
 304     def f_Rh(rh,h):
 305         return delta*rh*rh-v[h]*rh-sum([u[(h,i)]*eta[(h,i)] for i in N])
 306
 307     for h in V:
 308         if not ASYNC or  random() < taux_perte:
 309             rep = float(v[h])/(2*delta)
 310             Rhp[h] = 0 if rep < 0 else rep
 311         else :
 312             Rhp[h] = Rh[h]
 313
 314     xp={}
 315     def f_x(xhl,h,l):
 316         r =  delta*xhl*xhl + xhl*(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]))
 317         return r
 318
 319
 320     for h in V:
 321         for l in L:
 322             if not ASYNC or  random() < taux_perte:
 323                 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)
 324                 xp[(h,l)] = 0 if rep < 0 else rep
 325             else :
 326                 xp[(h,l)] = x[(h,l)]
 327
 328     valeurFonctionDualep = 0
 329     valeurFonctionDualep += sum([amplifieur*q[i]*q[i] for i in N])
 330     valeurFonctionDualep += sum([sum([delta*mt.pow(x[(h,l)],2) for l in L]) for h in V])
 331     valeurFonctionDualep += sum([delta*(Rh[h]**2) for h in V])
 332     valeurFonctionDualep += sum([sum([u[(h,i)]*(sum([ a[i][l]*x[(h,l)] for l in L])- eta[(h,i)]) for i in N]) for h in V])
 333     valeurFonctionDualep += sum([v[h]*(mt.log(float(sigma2)/D)/(gamma*mt.pow(Ps[h],float(2)/3)) - Rh[h]) for h in V])
 334     valeurFonctionDualep += sum([la[i]*(psi(Ps,i) +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]) -q[i]*Bi)  for i in N ])
 335     valeurFonctionDualep += sum([w[l]*sum([a[i][l]*q[i] for i in N]) for l in L])
 336
 337
 338     AfficheVariation(up,vp,lap,wp,thetap,etap,qp,Psp,Rhp,xp,valeurFonctionDualep)
 339
 340     arret = abs(valeurFonctionDuale-valeurFonctionDualep) < error
 341
 342     return (up,vp,lap,wp,thetap,etap,qp,Psp,Rhp,xp,valeurFonctionDualep,arret)
 343
 344
 345
 346
 347 """
 348 print "u",u
 349 print "v",v
 350 print "lambda",la
 351 print "w",w
 352 print "theta",theta
 353 print "eta", eta
 354 print "q",q
 355 print "Ps",Ps
 356 print "Rh",Rh
 357 print "x",x
 358
 359 """
 360
 361
 362
 363 def initialisation():
 364     global u, v, la, w, theta , q,  Ps, Rh,  eta, x, valeurFonctionDuale
 365
 366     theta = omega
 367
 368
 369     # TODO Ahmed : initialiser avec qi >0 et somme (ail.qi) = 0 Ei/i = Ti
 370     #q = [0 for i in N]
 371     q = {}
 372     for i in N :
 373         q[i] = 0.15 + random()*0.05
 374
 375
 376
 377
 378
 379     #TODO Ahmed doit donner les valeurs initiales pour Ps
 380     Ps= {}
 381     for h in V :
 382         Ps[h] =  0.2+random()*0.3
 383
 384
 385
 386
 387     #TODO Ahmed doit donner les valeurs initiales pour Rh
 388     Rh = {}
 389     for vi in V:
 390         Rh[vi] = 0.1 + random()*0.1
 391
 392
 393
 394     eta = {}
 395     for h in V :
 396         for i in N:
 397             eta[(h,i)]= etahi(h,i,Rh)
 398
 399
 400
 401     # Ok pour 0
 402     x={}
 403     for h in V :
 404         for l in L:
 405             x[(h,l)]=0#random()
 406
 407
 408
 409
 410     # initialisation des operateurs lagrangiens
 411     u={}
 412     for h in V :
 413         for i in N:
 414             u[(h,i)]= random()
 415
 416     v = {}
 417     for h in V :
 418         v[h] = Rh[h]
 419
 420     la = {}
 421     for i in N:
 422         la[i] = random()
 423
 424     w={}
 425     for l in L:
 426         w[l] =  random()
 427
 428
 429
 430
 431
 432 def __evalue_maj_theta():
 433     global u, v, la, w, theta , q,  Ps, Rh,  eta, x, valeurFonctionDuale
 434     nbexp = 1
 435     res = {}
 436     for l in range(1,11):
 437         m = 0
 438         def __maj_theta(k):
 439             om = omega/(mt.pow(k,float(l)/10))
 440             #om = omega/(mt.pow(k,0.5))
 441             #print "ici",k,om
 442             return om
 443         for _ in range(nbexp):
 444             initialisation()
 445             k = 1
 446             arret = False
 447             while k < 10000  and not arret :
 448                 (u,v,la,w,theta,eta,q,Ps,Rh,x,valeurFonctionDuale,ar)=maj(k,__maj_theta)
 449                 arret = ar
 450                 k+=1
 451             m += k
 452         res[l]= m/nbexp
 453         print l, res[l]
 454     print res
 455
 456
 457 __evalue_maj_theta()
 458 #ASYNC = True
 459 #__evalue_maj_theta()
 460
 461
 462 """
 463 initialisation()
 464 k = 1
 465 arret = False
 466 while k < 10000  and not arret :
 467     (u,v,la,w,theta,eta,q,Ps,Rh,x,valeurFonctionDuale,ar)=maj(k,maj_theta)
 468     arret = ar
 469     k+=1
 470 """
 471
 472
 473 """
 474     print "u",u
 475     print "w",w
 476     print "theta",theta
 477     print "eta", eta
 478     print "q",q
 479     print "v",v
 480     print "lambda",la
 481     print "Ps",Ps
 482     print "Rh",Rh
 483     print "x",x
 484 """
 485
 486
 487 """
 488     k +=1
 489
 490 """
 491
 492 """
 493 print "nb it",k
 494 print "u",u
 495 print "v",v
 496 print "lambda",la
 497 print "w",w
 498 print "theta",theta
 499 print "eta", eta
 500 print "q",q
 501 print "Ps",Ps
 502 print "Rh",Rh
 503 print "x",x
 504 print "L",valeurFonctionDuale
 505 """
 506
 507
 508 # relation ente les variables primaires et secondaires ?
 509