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 from copy import deepcopy
   9 error  = 0.1
  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 init = []
  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,mxg,idxexp):
 219     # initialisation des operateurs lagrangiens
 220     global u, v, la, w, theta , q,  Ps, Rh,  eta, x, valeurFonctionDuale
 221
 222     smax = 0
 223     #print "iteration",k
 224
 225     up = {}
 226     for h in V:
 227         for i in N:
 228             if not ASYNC or  random() < taux_succes:
 229                 s = eta[(h,i)]-sum([a[i][l]*x[(h,l)] for l in L])
 230                 if abs(s) > mxg :
 231                     print "ds calcul u",abs(s),idxexp
 232                     mxg = abs(s)
 233                 smax = max(smax,abs(s))
 234                 up[(h,i)] = u[(h,i)]-theta*s
 235             else :
 236                 up[(h,i)] = u[(h,i)]
 237
 238
 239
 240     vp = {}
 241     for h in V:
 242         if not ASYNC or  random() < taux_succes:
 243             s = Rh[h]- mt.log(float(sigma2)/D)/(gamma*mt.pow(Ps[h],float(1)/3))
 244             if abs(s) > mxg :
 245                 print "ds calcul v",abs(s),idxexp
 246                 mxg = abs(s)
 247             smax = max(smax,abs(s))
 248             vp[h] = max(0,v[h]-theta*s)
 249         else :
 250             vp[h] = v[h]
 251
 252
 253
 254     lap={}
 255     for i in N:
 256         if not ASYNC or  random() < taux_succes:
 257             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)
 258             if abs(s) > mxg :
 259                 print "ds calcul la",abs(s),idxexp,i
 260                 mxg = abs(s)
 261             smax = max(smax,abs(s))
 262             resla = la[i]-theta*s
 263             lap[i] = max(0,resla)
 264         else :
 265             lap[i] = la[i]
 266
 267
 268
 269
 270     wp={}
 271     for l in L:
 272         if not ASYNC or  random() < taux_succes:
 273             s = sum([a[i][l]*q[i] for i in N])
 274             if abs(s) > mxg :
 275                 print "ds calcul w",abs(s),idxexp
 276                 mxg = abs(s)
 277             smax = max(smax,abs(s))
 278             wp[l] = w[l] + theta*s
 279         else :
 280             wp[l] = w[l]
 281
 282
 283
 284     thetap = maj_theta(k)
 285
 286
 287     qp={}
 288     def f_q(qi,i):
 289         fa = sum([a[i][l]*w[l] for l in L]) - la[i]*Bi
 290         return qi*qi + qi*fa
 291
 292     for i in N:
 293         if not ASYNC or  random() < taux_succes:
 294             c = -float(sum([a[i][l]*w[l] for l in L]) - la[i]*Bi)/(2*amplifieur)
 295             rep = epsilon if c <= 0 else c
 296             qp[i] = rep
 297         else :
 298             qp[i] = q[i]
 299
 300
 301
 302
 303
 304     Psp={}
 305     #print "maj des des Psh"
 306     def f_Ps(psh,h):
 307         #print "ds f_ps",psh, v[h]* mt.log(float(sigma2)/D)/(gamma*((psh**2)**(float(1)/3))) +la[h]*psh
 308         return v[h]* mt.log(float(sigma2)/D)/(gamma*mt.pow(float(2)/3)) +la[h]*psh
 309     for h in V:
 310         if not ASYNC or  random() < taux_succes:
 311             lah = 0.05 if la[h] == 0 else  la[h]
 312             rep = (float(2*v[h]*mt.log(float(sigma2)/D))/mt.pow(3*gamma*lah,float(3)/5))
 313             Psp[h] = epsilon if rep <= 0 else rep
 314         else :
 315             Psp[h] = Ps[h]
 316
 317
 318
 319     etap={}
 320
 321     for h in V:
 322         for i in N :
 323             etap[(h,i)] = etahi(h,i,Rh)
 324
 325
 326
 327     Rhp={}
 328     def f_Rh(rh,h):
 329         return delta*rh*rh-v[h]*rh-sum([u[(h,i)]*eta[(h,i)] for i in N])
 330
 331     for h in V:
 332         if not ASYNC or  random() < taux_succes:
 333             rep = float(v[h])/(2*delta)
 334             Rhp[h] = 0 if rep < 0 else rep
 335         else :
 336             Rhp[h] = Rh[h]
 337
 338     xp={}
 339     def f_x(xhl,h,l):
 340         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]))
 341         return r
 342
 343
 344     for h in V:
 345         for l in L:
 346             if not ASYNC or  random() < taux_succes:
 347                 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)
 348                 xp[(h,l)] = 0 if rep < 0 else rep
 349             else :
 350                 xp[(h,l)] = x[(h,l)]
 351
 352
 353
 354
 355
 356     valeurFonctionDualep = 0
 357     valeurFonctionDualep += sum([amplifieur*q[i]*q[i] for i in N])
 358     valeurFonctionDualep += sum([sum([delta*(x[(h,l)]**2) for l in L]) for h in V])
 359     valeurFonctionDualep += sum([delta*(Rh[h]**2) for h in V])
 360     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])
 361     valeurFonctionDualep += sum([v[h]*(mt.log(float(sigma2)/D)/(gamma*mt.pow(Ps[h],float(2)/3)) - Rh[h]) for h in V])
 362     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 ])
 363     valeurFonctionDualep += sum([w[l]*sum([a[i][l]*q[i] for i in N]) for l in L])
 364
 365
 366     #AfficheVariation(up,vp,lap,wp,thetap,etap,qp,Psp,Rhp,xp,valeurFonctionDualep)
 367
 368     arret = abs(valeurFonctionDuale-valeurFonctionDualep) < error
 369
 370     return (up,vp,lap,wp,thetap,etap,qp,Psp,Rhp,xp,valeurFonctionDualep,arret,mxg,smax)
 371
 372
 373
 374
 375 """
 376 print "u",u
 377 print "v",v
 378 print "lambda",la
 379 print "w",w
 380 print "theta",theta
 381 print "eta", eta
 382 print "q",q
 383 print "Ps",Ps
 384 print "Rh",Rh
 385 print "x",x
 386
 387 """
 388
 389
 390
 391 def initialisation():
 392     global u, v, la, w, theta , q,  Ps, Rh,  eta, x,init
 393
 394     theta = omega
 395
 396     q = {}
 397     for i in N :
 398         q[i] = 0.15 + random()*0.05
 399
 400
 401     Ps= {}
 402     for h in V :
 403         Ps[h] =  0.2+random()*0.3
 404
 405     Rh = {}
 406     for vi in V:
 407         Rh[vi] = 0.1 + random()*0.1
 408
 409     eta = {}
 410     for h in V :
 411         for i in N:
 412             eta[(h,i)]= etahi(h,i,Rh)
 413
 414     x={}
 415     for h in V :
 416         for l in L:
 417             x[(h,l)]=0
 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     init = [deepcopy(q),deepcopy(Ps),deepcopy(Rh),deepcopy(eta),
 441             deepcopy(x),deepcopy(u),deepcopy(v),deepcopy(la),deepcopy(w)]
 442
 443
 444
 445 def __evalue_maj_theta__():
 446     global u, v, la, w, theta , q,  Ps, Rh,  eta, x, valeurFonctionDuale
 447     nbexp = 10
 448     res = {}
 449     m = []
 450     itermax = 100000
 451
 452     def __maj_theta(k):
 453         om = omega/(mt.pow(k,0.75))
 454         return om
 455     for idxexp in range(nbexp):
 456         mxg = 0
 457         initialisation()
 458         k = 1
 459         arret = False
 460         sm = 0
 461         while k < itermax  and not arret :
 462             (u,v,la,w,theta,eta,q,Ps,Rh,x,valeurFonctionDuale,ar,mxg,smax)=maj(k,__maj_theta,mxg,idxexp)
 463             errorq =  (max(q.values()) - min(q.values()))/min(q.values())
 464             arret = errorq <  error
 465             k+=1
 466             variation = "+" if smax > sm else "-"
 467             print variation,
 468             if k%100 ==0 :
 469                 print "k:",k,"erreur sur q", errorq, "et q:",q
 470                 print "maxg=", mxg
 471             if smax - sm  > 500:
 472                 print "variation trop grande"
 473                 print "init"
 474                 print init
 475                 exit
 476             sm = smax
 477
 478         if k == itermax:
 479             print "nbre d'iteration trop grand"
 480             print "init"
 481             print init
 482             exit
 483
 484         print "###############"
 485         print k
 486         print "###############"
 487         m += [k]
 488
 489     print (min(m),max(m),float(sum(m))/nbexp,m),m
 490
 491
 492
 493 def __une_seule_exp__(fichier_donees):
 494     global u, v, la, w, theta , q,  Ps, Rh,  eta, x, valeurFonctionDuale
 495     initialisation()
 496
 497
 498     fichier = open(fichier_donees, "r")
 499     instructions ={}
 500     for line in fichier:
 501          l = line.split("=")
 502          instructions[l[0]] = eval(l[1])
 503     u, v, la, w,  q,  Ps, Rh,  eta, x, = instructions['u'],    instructions['v'],    instructions['la'],    instructions['w'],    instructions['q'],    instructions['Ps'],    instructions['Rh'],    instructions['eta'],    instructions['x']
 504     nbexp = 1
 505     res = {}
 506     m = []
 507     itermax = 100000
 508
 509     def __maj_theta(k):
 510         om = omega/(mt.pow(k,0.75))
 511         return om
 512     for idxexp in range(nbexp):
 513         mxg = 0
 514         k = 1
 515         arret = False
 516         sm = 0
 517         while k < itermax  and not arret :
 518             (u,v,la,w,theta,eta,q,Ps,Rh,x,valeurFonctionDuale,ar,mxg,smax)=maj(k,__maj_theta,mxg,idxexp)
 519             errorq =  (max(q.values()) - min(q.values()))/min(q.values())
 520             arret = errorq <  error
 521             k+=1
 522             variation = "+" if smax > sm else "-"
 523             print variation,
 524             if k%100 ==0 :
 525                 print "k:",k,"erreur sur q", errorq, "et q:",q
 526                 print "maxg=", mxg
 527             if smax - sm  > 500:
 528                 print "variation trop grande"
 529                 print "init"
 530                 print init
 531                 exit
 532             sm = smax
 533
 534         if k == itermax:
 535             print "nbre d'iteration trop grand"
 536             print "init"
 537             print init
 538             exit
 539
 540         print "###############"
 541         print k
 542         print "###############"
 543         m += [k]
 544
 545     print (min(m),max(m),float(sum(m))/nbexp,m),m
 546
 547
 548
 549
 550 ASYNC = False
 551 __une_seule_exp__("config_initiale.py")
 552 #__evalue_maj_theta__()
 553 #ASYNC = True
 554 #taux_succes = 0.9
 555 #__evalue_maj_theta__()
 556
 557
 558
 559
 560
 561 print "u",u
 562 print "v",v
 563 print "lambda",la
 564 print "w",w
 565 print "theta",theta
 566 print "eta", eta
 567 print "q",q
 568 print "Ps",Ps
 569 print "Rh",Rh
 570 print "x",x
 571 print "L",valeurFonctionDuale
 572
 573
 574
 575 # relation ente les variables primaires et secondaires ?
 576