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 = 5
 448     res = {}
 449     m = []
 450     itermax = 10000
 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             sm = smax
 468             print variation,
 469             if k%100 ==0 :
 470                 print "k:",k,"erreur sur q", errorq, "et q:",q
 471                 print "maxg=", mxg
 472             if smax - sm  > 500:
 473                 print "variation trop grande"
 474                 print "init"
 475                 print init
 476                 exit
 477         if k == itermax:
 478             print "nbre d'iteration trop grand"
 479             print "init"
 480             print init
 481             exit
 482
 483         print "###############"
 484         print k
 485         print "###############"
 486         m += [k]
 487
 488     print (min(m),max(m),float(sum(m))/nbexp,m),m
 489
 490
 491
 492
 493
 494 ASYNC = False
 495 __evalue_maj_theta__()
 496 #ASYNC = True
 497 #taux_succes = 0.9
 498 #__evalue_maj_theta__()
 499
 500
 501 """
 502 initialisation()
 503 k = 1
 504 arret = False
 505 while k < 10000  and not arret :
 506     (u,v,la,w,theta,eta,q,Ps,Rh,x,valeurFonctionDuale,ar)=maj(k,maj_theta)
 507     arret = ar
 508     k+=1
 509     errorq =  abs(min(q.values())-max(q.values()))
 510     print "errorq",errorq
 511     arret = errorq <  error
 512 """
 513
 514 """
 515     print "u",u
 516     print "w",w
 517     print "theta",theta
 518     print "eta", eta
 519     print "q",q
 520     print "v",v
 521     print "lambda",la
 522     print "Ps",Ps
 523     print "Rh",Rh
 524     print "x",x
 525 """
 526
 527
 528 """
 529     k +=1
 530
 531 """
 532
 533
 534
 535 print "u",u
 536 print "v",v
 537 print "lambda",la
 538 print "w",w
 539 print "theta",theta
 540 print "eta", eta
 541 print "q",q
 542 print "Ps",Ps
 543 print "Rh",Rh
 544 print "x",x
 545 print "L",valeurFonctionDuale
 546
 547
 548
 549 # relation ente les variables primaires et secondaires ?
 550