exp_controle_asynchrone/simulMWSN.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  = 0.01
  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
  21
  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 #print "L",L
  78 #print "d",d
  79
  80
  81 def ail(i,l):
  82     assert  l in L, " pb de dimennsion de l: "+str(l)+ " "+ str(L)
  83     r = 0
  84     (o,e) = G.edges()[l]
  85     if o == i :
  86         r = 1
  87     elif e == i :
  88         r = -1
  89     return r
  90
  91
  92
  93 # Constantes
  94 a = [[ail(i,l) for l in L ] for i in xrange(n)]
  95 aplus  = [[1 if  ail(i,l) > 0 else 0  for l in L ] for i in xrange(n)]
  96 amoins  = [[1 if  ail(i,l) < 0 else 0  for l in L ] for i in xrange(n)]
  97
  98
  99
 100
 101
 102 alpha = 0.5
 103 beta = 1.3E-8
 104 gamma = 55.54
 105 delta = 0.2
 106 amplifieur = 1
 107 sigma2 = 3500
 108 Bi = 5
 109 omega = 0.15
 110 D = 100
 111 path_loss_exp = 4
 112 cr = 0.5
 113 cs = [alpha + beta*(di**path_loss_exp) for di in d]
 114 #print "cs",cs
 115
 116
 117 #TODO
 118 def etahi(h,i,R):
 119     ret = 0
 120     if i == h :
 121         ret = R[h]
 122     elif i == sink  :
 123         ret = -R[h]
 124     return ret
 125
 126
 127 def psi(Ps,i):
 128     if i not in Ps:
 129         return 0
 130     else :
 131         return Ps[i]
 132
 133
 134
 135 def cmpPair(x1,x2):
 136     return cmp(x1[1],x2[1])
 137
 138 def aminp(l):
 139     l.sort(cmp=cmpPair)
 140     return l[0][0]
 141
 142
 143
 144 def distance(d1,d2):
 145     return mt.sqrt(sum([(d1[t]-d2[t])**2 for t in d1]))
 146
 147
 148 def AfficheVariation (up,vp,lap,wp,thetap,etap,qp,Psp,Rhp,xp,valeurFonctionDualep):
 149     global u, v, la, w, theta , q, Ps, Rh, eta, x,valeurFonctionDuale
 150
 151     print "du=",distance(u,up),
 152     print "dv=",distance(v,vp),
 153     print "dlambda=",distance(la,lap),
 154     print "dw=",distance(w,wp),
 155     print "dtheta=",abs(theta-thetap),
 156     print "deta=",distance(eta,etap),
 157     print "dq=",distance(q,qp),
 158     print "dPs=",distance(Ps,Psp),
 159     print "dRh=",distance(Rh,Rhp),
 160     print "dx=",distance(x,xp),
 161     print "dL=",abs(valeurFonctionDuale-valeurFonctionDualep),"\n"
 162
 163
 164 valeurFonctionDuale = 0
 165
 166 def entre0et1(x):
 167     r = x if (x >0 and x <= 1) else -1
 168     #print "ds entre0et1 x et r", x,r
 169     return r
 170
 171
 172 def xpos(x):
 173     r = x if x >0 else -1
 174     #print "ds xpos x et r", x,r
 175     return r
 176
 177 def xposounul(x):
 178     return x if x >= 0 else -1
 179
 180 #f_q,q[i],POS,[i]
 181 def armin(f,xini,xr,args):
 182     #xpos = lambda x : x if x > 0 else -1
 183     r = 0
 184     if xr == POS :
 185         #print "strictement pos"
 186         #print "parametres passes a cobyla",xini,xpos,args,"consargs=(),rhoend=1E-5,iprint=0,maxfun=1000"
 187         r= opt.fmin_cobyla(f,xini,cons=[xpos],args=args,consargs=(),rhoend=1E-3,iprint=0,maxfun=nbiter)
 188         #print "le min str pos est",r
 189         #print "le min pos   est", r,"avec xini",xini
 190     elif xr == POS_NUL :
 191         #print "pos ou nul"
 192         r = opt.fmin_cobyla(f,xini,[xposounul],args,consargs=(),rhoend=1E-3,iprint=0,maxfun=nbiter)
 193         # print "le min pos est", r
 194         #print "le min pos null  est", r,"avec xini",xini
 195     elif xr == POSINF1:
 196         r = opt.fmin_cobyla(f,xini,[entre0et1],args,consargs=(),rhoend=1E-3,iprint=0,maxfun=nbiter)
 197         #print "le min pos inf 1 est", r,"avec xini",xini
 198
 199
 200     return r
 201
 202
 203
 204
 205
 206
 207 def maj_theta(k):
 208     return omega/(mt.pow(k,0.5))
 209
 210
 211
 212 def maj_theta(k):
 213     return omega/mt.sqrt(k)
 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     vp = {}
 236     for h in V:
 237         #print "vp",gamma*mt.pow(Ps[h],float(1)/3)
 238         #vp[h] = max(0,v[h]-theta*(Rh[h]- mt.log(float(sigma2)/D)/(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
 246     lap={}
 247     for i in N:
 248         if not ASYNC or  random() < taux_perte:
 249             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)))
 250             lap[i] = max(0,resla)
 251         else :
 252             lap[i] = la[i]
 253
 254
 255
 256
 257     wp={}
 258     for l in L:
 259         if not ASYNC or  random() < taux_perte:
 260             wp[l] = w[l] + theta*(sum([a[i][l]*q[i] for i in N]))
 261         else :
 262             wp[l] = w[l]
 263
 264
 265
 266     thetap = maj_theta(k)
 267
 268     #print "k",k,"maj theta", thetap
 269
 270
 271     qp={}
 272     def f_q(qi,i):
 273         fa = sum([a[i][l]*w[l] for l in L]) - la[i]*Bi
 274         return qi*qi + qi*fa
 275
 276     for i in N:
 277         if not ASYNC or  random() < taux_perte:
 278             c = -float(sum([a[i][l]*w[l] for l in L]) - la[i]*Bi)/(2*amplifieur)
 279             rep = epsilon if c <= 0 else c
 280             qp[i] = rep
 281         else :
 282             qp[i] = q[i]
 283
 284
 285
 286
 287
 288     Psp={}
 289     #print "maj des des Psh"
 290     def f_Ps(psh,h):
 291         #print "ds f_ps",psh, v[h]* mt.log(float(sigma2)/D)/(gamma*((psh**2)**(float(1)/3))) +la[h]*psh
 292         return v[h]* mt.log(float(sigma2)/D)/(gamma*mt.pow(float(2)/3)) +la[h]*psh
 293     for h in V:
 294         if not ASYNC or  random() < taux_perte:
 295             lah = 0.05 if la[h] == 0 else  la[h]
 296             rep = (float(2*v[h]*mt.log(float(sigma2)/D))/mt.pow(3*gamma*lah,float(3)/5))
 297             Psp[h] = epsilon if rep <= 0 else rep
 298         else :
 299             Psp[h] = Ps[h]
 300
 301
 302
 303     etap={}
 304
 305     for h in V:
 306         for i in N :
 307             etap[(h,i)] = etahi(h,i,Rh)
 308
 309
 310
 311     Rhp={}
 312     def f_Rh(rh,h):
 313         return delta*rh*rh-v[h]*rh-sum([u[(h,i)]*eta[(h,i)] for i in N])
 314
 315     for h in V:
 316         if not ASYNC or  random() < taux_perte:
 317             rep = float(v[h])/(2*delta)
 318             Rhp[h] = 0 if rep < 0 else rep
 319         else :
 320             Rhp[h] = Rh[h]
 321
 322     xp={}
 323     def f_x(xhl,h,l):
 324         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]))
 325         return r
 326
 327
 328     for h in V:
 329         for l in L:
 330             if not ASYNC or  random() < taux_perte:
 331                 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)
 332                 xp[(h,l)] = 0 if rep < 0 else rep
 333             else :
 334                 xp[(h,l)] = x[(h,l)]
 335
 336
 337
 338
 339
 340     valeurFonctionDualep = 0
 341     valeurFonctionDualep += sum([amplifieur*q[i]*q[i] for i in N])
 342     valeurFonctionDualep += sum([sum([delta*(x[(h,l)]**2) for l in L]) for h in V])
 343     valeurFonctionDualep += sum([delta*(Rh[h]**2) for h in V])
 344     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])
 345     valeurFonctionDualep += sum([v[h]*(mt.log(float(sigma2)/D)/(gamma*mt.pow(Ps[h],float(2)/3)) - Rh[h]) for h in V])
 346     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 ])
 347     valeurFonctionDualep += sum([w[l]*sum([a[i][l]*q[i] for i in N]) for l in L])
 348
 349
 350     #AfficheVariation(up,vp,lap,wp,thetap,etap,qp,Psp,Rhp,xp,valeurFonctionDualep)
 351
 352     arret = abs(valeurFonctionDuale-valeurFonctionDualep) < error
 353
 354     return (up,vp,lap,wp,thetap,etap,qp,Psp,Rhp,xp,valeurFonctionDualep,arret)
 355
 356
 357
 358
 359 """
 360 print "u",u
 361 print "v",v
 362 print "lambda",la
 363 print "w",w
 364 print "theta",theta
 365 print "eta", eta
 366 print "q",q
 367 print "Ps",Ps
 368 print "Rh",Rh
 369 print "x",x
 370
 371 """
 372
 373
 374
 375 def initialisation():
 376     global u, v, la, w, theta , q,  Ps, Rh,  eta, x, valeurFonctionDuale
 377
 378     theta = omega
 379
 380
 381     # TODO Ahmed : initialiser avec qi >0 et somme (ail.qi) = 0 Ei/i = Ti
 382     #q = [0 for i in N]
 383     q = {}
 384     for i in N :
 385         q[i] = 0.15 + random()*0.05
 386
 387
 388
 389
 390
 391     #TODO Ahmed doit donner les valeurs initiales pour Ps
 392     Ps= {}
 393     for h in V :
 394         Ps[h] =  0.2+random()*0.3
 395
 396
 397
 398
 399     #TODO Ahmed doit donner les valeurs initiales pour Rh
 400     Rh = {}
 401     for vi in V:
 402         Rh[vi] = 0.1 + random()*0.1
 403
 404
 405
 406     eta = {}
 407     for h in V :
 408         for i in N:
 409             eta[(h,i)]= etahi(h,i,Rh)
 410
 411
 412
 413     # Ok pour 0
 414     x={}
 415     for h in V :
 416         for l in L:
 417             x[(h,l)]=0#random()
 418
 419
 420
 421
 422     # initialisation des operateurs lagrangiens
 423     u={}
 424     for h in V :
 425         for i in N:
 426             u[(h,i)]= random()
 427
 428     v = {}
 429     for h in V :
 430         v[h] = Rh[h]
 431
 432     la = {}
 433     for i in N:
 434         la[i] = random()
 435
 436     w={}
 437     for l in L:
 438         w[l] =  random()
 439
 440
 441
 442 def __evalue_maj_theta__():
 443     global u, v, la, w, theta , q,  Ps, Rh,  eta, x, valeurFonctionDuale
 444     nbexp = 50
 445     res = {}
 446     m = []
 447     def __maj_theta(k):
 448         om = omega/(mt.pow(k,0.75))
 449         return om
 450     for idxexp in range(nbexp):
 451
 452         initialisation()
 453         k = 1
 454         arret = False
 455         while k < 1000000  and not arret :
 456             (u,v,la,w,theta,eta,q,Ps,Rh,x,valeurFonctionDuale,ar)=maj(k,__maj_theta)
 457             errorq =  abs(min(q.values())-max(q.values()))
 458             arret = errorq <  error
 459             k+=1
 460             if k%5000 ==0 :
 461                 print "k:",k,"erreur sur q", errorq, "et q:",q
 462         print "###############"
 463         print k
 464         print "###############"
 465         m += [k]
 466     print (min(m),max(m),float(sum(m))/nbexp,m),m
 467
 468
 469
 470
 471
 472 #ASYNC = False
 473 #__evalue_maj_theta()
 474 ASYNC = True
 475 taux_perte = 0.9
 476 __evalue_maj_theta__()
 477
 478
 479 """
 480 initialisation()
 481 k = 1
 482 arret = False
 483 while k < 10000  and not arret :
 484     (u,v,la,w,theta,eta,q,Ps,Rh,x,valeurFonctionDuale,ar)=maj(k,maj_theta)
 485     arret = ar
 486     k+=1
 487     errorq =  abs(min(q.values())-max(q.values()))
 488     print "errorq",errorq
 489     arret = errorq <  error
 490 """
 491
 492 """
 493     print "u",u
 494     print "w",w
 495     print "theta",theta
 496     print "eta", eta
 497     print "q",q
 498     print "v",v
 499     print "lambda",la
 500     print "Ps",Ps
 501     print "Rh",Rh
 502     print "x",x
 503 """
 504
 505
 506 """
 507     k +=1
 508
 509 """
 510
 511
 512
 513 print "u",u
 514 print "v",v
 515 print "lambda",la
 516 print "w",w
 517 print "theta",theta
 518 print "eta", eta
 519 print "q",q
 520 print "Ps",Ps
 521 print "Rh",Rh
 522 print "x",x
 523 print "L",valeurFonctionDuale
 524
 525
 526
 527 # relation ente les variables primaires et secondaires ?
 528