full simu with diag

algo.py

@@ -28,7 +28,10 @@ def __init__(self, DSSObj, env, P, Q, slot_len_in_min=10, start_hr=0, mode='chan
 
         #self.max_rate = 1+60.0*np.array(env['demand']) # In kWmin
         self.max_rate = 8.0*np.ones(env['evNumber']) # In kW
+        self.max_rate_scaler = 1.0
 
+        #self.remaining_demand /= self.max_rate_scaler
+
         self.arrival = np.round(np.array(self.env['evArrival'])/slot_len_in_min)
         self.duration = np.round(np.array(self.env['evDuration'])/slot_len_in_min)
         self.claimed_duration = np.round(np.array(self.env['evClaimedDuration'])/slot_len_in_min)
@@ -105,7 +108,7 @@ def get_trans_load(self, ev_power, P, Q): # In kVA
 
     def get_available(self, trans_loads):
 
-        return self.trans_accu*np.maximum(0.0, np.array(self.env['transRating']) - trans_loads)
+        return self.trans_accu*np.maximum(0.0, np.array(self.env['transRating']) - trans_loads)/self.max_rate_scaler
 
     def get_UB(self, connected):
         #temp = np.minimum(self.remaining_demand/self.slot_len_in_min, self.max_rate)

central_algo.py

@@ -40,10 +40,11 @@ def update(self, P, Q):
             #x = lb.solve(self.get_driver_type(connected), self.get_claimed(connected), self.get_over_time(connected), self.get_UB(connected), A, T)
             #w = util.f(self.get_driver_type(connected)) * util.g(self.get_laxity(connected, scale=0.0))
             w = util.w(self.get_discrepancy(connected), self.get_laxity(connected, scale=0.0))
+            #print(w)
             #laxity = self.get_laxity(connected)
             #w = w * (144-laxity)
 
-            x = lb.solve(w=w, UB=self.get_UB(connected), A=A, T=T)
+            x = lb.solve(w=w, UB=self.get_UB(connected), A=A, T=T, scale=self.max_rate_scaler)
             #print(x)
             #x = primal.solve(np.zeros(len(connected)), 10*np.ones(len(connected)), A, T)
             #x = primal.solve(np.zeros(len(connected)), self.get_UB(connected)-LB+util.tol, A, T)

conv_ana_3.py

@@ -6,7 +6,13 @@
 
 from tqdm import tqdm
 from matplotlib import pyplot as plt
+'''
+import sys
 
+i = 1
+print(str(i).zfill(5))
+sys.exit()
+'''
 result_path = 'result/test_500.txt'
 base_load_path = 'base_load/10_min/'
 env_path = 'env/test_500.txt'
@@ -18,13 +24,13 @@
 result = util.load_dict(result_path)
 env = util.load_dict(env_path)
 
-slot = 102
+slot = 109
 rho = 1000.0
 #print(result['central'][slot]['x'])
 connected = result['central'][slot]['connected']
 print('connected')
 print(len(connected))
-factor = 0.3
+factor = 0.1
 
 n_slot_per_hr = 6
 h = slot//n_slot_per_hr
@@ -42,7 +48,7 @@ def get_driver_type():
     return np.array(driver_type)
 
 #w = util.f(get_driver_type())
-w = result['central'][slot]['w']
+w = np.array(result['central'][slot]['w'])+10000
 
 def get_trans_load(ev_power): # In kVA
     DSSCircuit = RunPF.runPF(DSSObj, P[:, h%n_slot_per_hr], Q[:, h%n_slot_per_hr], env['evNodeNumber'], ev_power)
@@ -55,6 +61,7 @@ def get_trans_load(ev_power): # In kVA
     trans_loads = [np.sqrt(trans_loads[i]**2+trans_loads[i+1]**2) for i in range(0,len(trans_loads),2)]
     return(np.array(trans_loads))
 
+    #        connected_ = n
 def get_UB():
     return 8.0*np.ones(len(connected))
 
@@ -98,7 +105,8 @@ def get_TAU(whole=0):
 
 def get_load_nabla(T,A,U,LB,mu):
 
-    x = np.minimum(np.maximum(LB, w/util.non_zero(np.dot(T.T, mu**2))), get_UB())
+    #x = np.minimum(np.maximum(LB, w/util.non_zero(np.dot(T.T, mu**2))), get_UB())
+    x = np.minimum(np.maximum(LB, w/util.non_zero(np.dot(T.T, mu))), get_UB())
 
     ev_power = np.zeros(env['evNumber'])
     for j in range(0, len(connected)):
@@ -111,7 +119,8 @@ def get_load_nabla(T,A,U,LB,mu):
 
     rating_load = np.array([nabla[e] for e in U])
 
-    nabla = 2 * mu * rating_load
+    #nabla = 2 * mu * rating_load
+    nabla = rating_load
 
     return (load, nabla, x, rating_load)
 
@@ -125,21 +134,25 @@ def get_mu(mu_k, load_k, load_k_1, nabla_k, nabla_k_1, gamma=None):
 '''
 def get_mu(mu_k, mu_k_1, load_k, load_k_1, nabla_k, rating_load, gamma=None):
     #rating_load = nabla_k / (mu_k+tol)
-    hessian = rating_load - 2 * mu_k * ((load_k-load_k_1)/util.non_zero((mu_k-mu_k_1+tol)))
+    #hessian = rating_load - 2 * mu_k * ((load_k-load_k_1)/util.non_zero((mu_k-mu_k_1+tol)))
+    hessian = -1.0*((load_k-load_k_1) / util.non_zero(mu_k-mu_k_1))
     #hessian = np.array([util.tol if e <= util.tol else e for e in hessian])
     hessian = factor / util.non_zero(hessian) 
     if gamma==None: 
-        mu = mu_k - hessian * nabla_k
+        #mu = mu_k - hessian * nabla_k
+        #mu = np.maximum(0.0, mu_k - hessian * nabla_k)
+        mu = np.abs(mu_k - hessian * nabla_k)
     else:
         #mu = rating_load - 2 * mu_k * ((load_k-load_k_1)/(mu_k-mu_K_1))
-        mu = mu_k - gamma * hessian * nabla_k
+        #mu = mu_k - gamma * hessian * nabla_k
+        mu = np.maximum(0.0, mu_k - gamma * hessian * nabla_k)
     return mu
 
 T, A, U = get_TAU()
-scale = 1e-4
+scale = 1e-3
 legend = []
 
-for name in ['gpa']:
+for name in ['diag']:
     gammas = []
     iters = []
     for gamma in tqdm(range(5, 100)):
@@ -157,7 +170,7 @@ def get_mu(mu_k, mu_k_1, load_k, load_k_1, nabla_k, rating_load, gamma=None):
 
         if name == 'diag':
             ##################################################################
-            mu_k_1 = 0.0*np.ones(len(A))
+            mu_k_1 = np.zeros(len(A))
             load_k_1 = np.zeros(len(A))
 
             mu_k = np.ones(len(A))
@@ -198,7 +211,7 @@ def get_mu(mu_k, mu_k_1, load_k, load_k_1, nabla_k, rating_load, gamma=None):
                 mu_k_1 = np.copy(mu_k)
 
                 (load_k, nabla_k, x, rating_load) = get_load_nabla(T,A,U,LB,mu)
-
+                #print(mu_k)
                 mu_k = np.copy(mu)
 
                 #lamda = lamda_k - ( ( load_k-load_k_1 )/( nabla_k-nabla_k_1 ) ) * nabla_k
@@ -264,15 +277,15 @@ def get_mu(mu_k, mu_k_1, load_k, load_k_1, nabla_k, rating_load, gamma=None):
             sub[1] = temp[0]+temp[1]+temp[2]
             
             #print(abs(sub[1]-sub[0])/sub[0])
-            if abs(sub[1]-sub[0]) <= 0.05*sub[0]:
+            if abs(sub[1]-sub[0]) <= 0.20*sub[0]:
                 break
             
             '''
 
             c = sum(result['central'][slot]['ev_power'])
             d = sum(ev_power)
 
-            if abs(d-c) <= 0.05*c:
+            if abs(d-c) <= 0.01*c:
                 break 
 
             #if np.allclose(x, result['central'][slot]['x'], atol=0.0, rtol=0.05)==True:

diag_algo.py

@@ -7,7 +7,8 @@
 class diag_algo(algo):
     def get_load_nabla(self, T,A,U,LB,UB,connected,P,Q,w,mu):
 
-        x = np.minimum(np.maximum(LB, w/util.non_zero(np.dot(T.T, mu**2))), UB)
+        #x = np.minimum(np.maximum(LB, w/util.non_zero(np.dot(T.T, mu**2))), UB)
+        x = np.minimum(np.maximum(LB, w/util.non_zero(np.dot(T.T, mu))), UB)
 
         ev_power = np.zeros(self.env['evNumber'])
         for j in range(0, len(connected)):
@@ -20,17 +21,21 @@ def get_load_nabla(self, T,A,U,LB,UB,connected,P,Q,w,mu):
 
         rating_load = np.array([nabla[e] for e in U])
 
-        nabla = 2 * mu * rating_load
+        #nabla = 2 * mu * rating_load
+        nabla = rating_load
 
         return (load, nabla, x, rating_load)
 
     def get_mu(self, mu_k, mu_k_1, load_k, load_k_1, nabla_k, rating_load, gamma=None):
         #rating_load = nabla_k / (mu_k+tol)
-        hessian = np.absolute(rating_load - 2 * mu_k * ((load_k-load_k_1)/util.non_zero((mu_k-mu_k_1+0.5))))
+        #hessian = np.absolute(rating_load - 2 * mu_k * ((load_k-load_k_1)/util.non_zero((mu_k-mu_k_1+0.5))))
+        hessian =  0.1+np.absolute( ( (load_k-load_k_1)/util.non_zero((mu_k-mu_k_1)) ) )
+        #print(hessian)
         #hessian = np.array([util.tol if e <= util.tol else e for e in hessian])
-        hessian = 0.3 / util.non_zero(hessian) 
+        hessian = 0.1 / util.non_zero(hessian) 
         if gamma==None: 
-            mu = mu_k - hessian * nabla_k
+            #mu = mu_k - hessian * nabla_k
+            mu = np.maximum(0,  mu_k - hessian * nabla_k )
         else:
             #mu = rating_load - 2 * mu_k * ((load_k-load_k_1)/(mu_k-mu_K_1))
             mu = mu_k - gamma * hessian * nabla_k
@@ -72,10 +77,10 @@ def update(self, P, Q, central={}):
 
             gamma = 0.008
 
-            mu_k_1 = 0.0*np.ones(len(A))
+            mu_k_1 = 1000*np.ones(len(A))
             load_k_1 = np.zeros(len(A))
 
-            mu_k = 1.0*np.ones(len(A))
+            mu_k = 999*np.ones(len(A))
 
             (load_k, nabla_k, _, rating_load) = self.get_load_nabla(T,A,U,LB,UB,connected,P,Q,w,mu_k)
             #print('gamma')
@@ -94,6 +99,7 @@ def update(self, P, Q, central={}):
                 mu_k_1 = np.copy(mu_k)
 
                 (load_k, nabla_k, x, rating_load) = self.get_load_nabla(T,A,U,LB,UB,connected,P,Q,w,mu)
+                x *= self.max_rate_scaler
 
                 mu_k = np.copy(mu)
 

gpa_algo.py

@@ -39,7 +39,7 @@ def update(self, P, Q, central={}):
             #gamma = (2.0*self.params['step_factor'])/(m*L*S+util.tol)
             #gamma = 0.00060483158055174284
 
-            gamma = 0.0015
+            gamma = 0.0000001
             #print('gamma')
             #print(gamma)
             #lamda = np.zeros(len(A))

lower_bound.py

@@ -2,14 +2,14 @@
 import cvxpy as cp
 import utility as util
 
-def solve(w, UB, A, T):
+def solve(w, UB, A, T, scale=1.0):
     #w = util.f(driver_type)
 
     z = cp.Variable(len(UB))
 
     #print(A)
     #print(w)
-    obj = cp.Maximize(sum(np.diag(w) * cp.log(z + util.tol)))
+    obj = cp.Maximize(sum(np.diag(w) * cp.log(100*z + util.tol)))
     #print(T)
     #print(A)
     constraints = [0 <= z, z <= UB, T * z <= A]
@@ -19,4 +19,4 @@ def solve(w, UB, A, T):
     # Solve with MOSEK.
     prob.solve(solver=cp.MOSEK, verbose=False)
     #print(z.value)
-    return z.value 
+    return z.value * scale

plot.py

@@ -323,7 +323,7 @@ def fig_conv_ana(result_path):
 
     #fig_soc_vs_time(simu_params['save_path'], (env['evDriverType']), algo='central')
     #fig_trans_load_vs_time(simu_params['save_path'], trans=0, env=env)
-    fig_trans_load_subplot(simu_params['save_path'], trans_list=[2,1,0], env=env)
-    #fig_compare(simu_params['save_path'], 1, 143, env)
+    #fig_trans_load_subplot(simu_params['save_path'], trans_list=[2,1,0], env=env)
+    fig_compare(simu_params['save_path'], 1, 143, env)
     #fig_conv_ana('result/meta_large.txt')