Secure MPC: Pressure Control Algorithm

111 days ago by ktjell12@student.aau.dk

##################################################################################################### # This cell implementes most of the secure MPC funtions, used in the master's thesis # # by Katrine Tjell. # # To name a few: # # - Shamir's Secret sharing scheme # # - Multiplication protocol # # - Beavers triplet creation # # - Integer division # # - Comparison protocol # ##################################################################################################### import numpy as np import matplotlib.pylab as plt import csv import sys def basispoly(F, n): r = [] C = range(1,n+1) for i in range(1,n+1): c = [k for k in C if k != i] p = 1 for j in range(n-1): p *= -F(c[j]) / (F(i)-F(c[j])) r.append(F(p)); return r def secretsharing(F, x, t, n): shares = [] c = [F.random_element() for i in range(t)] for i in range(1, n+1): s = x for j in range(1,t+1): s = F(s) + F(c[j-1]) * F(i)**F(j) shares.append(s) return shares def triplet(F, n, t): r = basispoly(F, n) matrixA = np.zeros((n,n)) matrixB = np.zeros((n,n)) for i in range(n): a = F.random_element() b = F.random_element() matrixA[:,i] = secretsharing(F, a, t, n) matrixB[:,i] = secretsharing(F, b, t, n) sharesA = np.sum(matrixA, 1) sharesA = [F(i) for i in sharesA] sharesB = np.sum(matrixB, 1) sharesB = [F(i) for i in sharesB] sharesC = [F(i*j) for (i,j) in zip(sharesA, sharesB)] return [sharesA, sharesB, sharesC] def dot(F, x, y): res = 0 for i in range(len(x)): res += F(x[i] * y[i]) return res #Should and could be implemented securely. def randomBit(F, n, t): #basis = basispoly(F,n) #rsq = 0 #while rsq == 0: #If r**2 = 0, the parties retry # r_shares = triplet(F,n,t)[0] #The parties compute an unknown random number, r, and computes r**2 # salpha, sbeta, sgamma = triplet(F, n, t) # d = dot(F, basis, [F(i-j) for (i,j) in zip(r_shares, salpha)]) # e = dot(F, basis, [F(i-j) for (i,j) in zip(r_shares, sbeta)]) # rsq_shares = [F(d*e + d*i + e*j + k) for (i,j,k) in zip(sbeta, salpha, sgamma)] # rsq = dot(F, rsq_shares, basis) #r**2 is revealed #r_prime = sqrt(rsq) #print 'rprime', r_prime r = np.random.choice(2,1)[0] return list(secretsharing(F,r,t,n))#list(F(1/2) * (1/r_prime * np.array(r_shares) + 1)) def randomBitNumber(F, n, t, l): basis = basispoly(F,n) sharings = [] for i in range(l): sharings.append(randomBit(F,n,t)) sharings_reverse = sharings[::-1] r_shares = [] for i in range(n): s = 0 for j in range(l): s += 2**j * dot(F, basis, sharings_reverse[j]) r_shares.append(s) return sharings, r_shares def mult(sx1,sx2,): salpha, sbeta, sgamma = triplet(F, n, t) sx1 = np.array(sx1) sx2 = np.array(sx2) r = basispoly(F,len(sx1)) d = dot(F, r, [F(i-j) for (i,j) in zip(sx1, salpha)]) e = dot(F, r, [F(i-j) for (i,j) in zip(sx2, sbeta)]) return np.array([F(d*e + d*i + e*j + k) for (i,j,k) in zip(sbeta, salpha, sgamma)]) def xor(a1,b1): d = np.array(a1) - np.array(b1) return mult(d,d) def OR(a): p = 1 - a[0] for i in range(1,len(a)): p = mult(p, (1 - a[i])) return 1 - p def lessthan(a,b): e = [] for i in range(len(a)): e.append(xor(a[i],b[i])) f = [] for i in range(len(e)): f.append(OR(e[0:i+1])) g = [f[0]] for i in range(0,len(f)-1): g.append(f[i+1] - f[i]) h = [] for i in range(len(g)): h.append(mult(g[i], b[i])) y = [] for i in range(len(h[0])): s = 0 for j in range(len(h)): s += h[j][i] y.append(s) return y def maxV(x): h = x[0] for i in range(1,len(x)): c = lessthan(h,x[i]) k = [1 - j for j in c] t1 = [mult(j,k) for j in h] t2 = [mult(c,j) for j in x[i]] h = [i+j for i,j in zip(t1,t2)] return h def berlekampWelch(F, x, y, t): b = np.multiply(y, x**t) v = vector(F, tuple(b)) A = [] for i in range(t): A.append(list(np.multiply(y, x**i))) for i in range(2*t+1): A.append(list(-x**i)) print('Shape of A: {}'.format(np.shape(A))) M = Matrix(F, A) M = M.transpose() print(M.echelon_form()) c = M.solve_right(v) print('x: {}'.format(c)) e = c[:t] h = c[t:] return e,h def bit2num(x): t = [] l = len(x) n = len(x[0]) for i in range(l): t.append([2**(l-1-i)*j for j in x[i]]) res = np.zeros(n, dtype = int) for i in range(l): res = np.array(res) + np.array(t[i], dtype = int) return res ###############################MAIN ########################################## ##### CONSTANTS ################################ m = 1125899839733759 # Prime number, large enough that the probability of chossing a certain number is negeligable F = GF(m) # Finite field with m elements n = 4 # Number of parties t = 1 # Degree of polynomial l = 8 # bit length of values x = [20,2] # Secret input values basis = basispoly(F,n) # Lagrange basispolynomial coefficients 
       
########################################################################################### # This cells computes SECURELY the pressure control algortihm # # - using the functions from cell 1 # #(thus cell 1 must run before attempting to run this cell). # ########################################################################################### import numpy as np import matplotlib.pylab as plt np.random.seed(1) k = 10 #Period length r = vector(k*[4]) #Pressure requirement T = 5 l = 4 c = [[0,0,0,0,0,3,3,3,3,3]] c1 = [[1,1,1,1,1,2,2,2,2,2]] cons = np.array([np.repeat(c, T, axis = 0).flatten(), np.repeat(c1, T, axis = 0).flatten()]) E = np.zeros((T, k)) Es = [] Y1 = np.zeros((T,k)) Y2 = np.zeros((T,k)) kappa = vector(k*[4]) kappaS = [] for i in range(k): kappaS.append(secretsharing(F,kappa[i],t,n)) for i in range(T): kappa = [] for j in kappaS: kappa.append(dot(F,basis,j)) Y = np.array([kappa - j for j in cons[:,i*k:k*(i+1)]]) Y1[i] = Y[0,:] Y2[i] = Y[1,:] error = np.array([r - j for j in Y]) s1 = [] s2 = [] for j in range(k): if error[0][j] == -1: error[0][j] = 0 if error[1][j] == -1: error[1][j] = 0 sb1 = bin(error[0][j])[2:].rjust(l,'0') sb2 = bin(error[1][j])[2:].rjust(l,'0') t1 = [] t2 = [] for q in range(l): t1.append(secretsharing(F,sb1[q],t,n)) t2.append(secretsharing(F,sb2[q],t,n)) s1.append(t1) s2.append(t2) em = [] for j in range(k): em.append(maxV([s1[j], s2[j]])) e_ks = [bit2num(j) for j in em] t1 = [] for j in kappaS: t1.append(dot(F,basis,j)) Es.append(t1) kappaS = np.array(kappaS) + np.array(e_ks) Es = np.array(Es) g = Graphics() g += list_plot(cons[0], plotjoined = True, legend_label = 'Consump1') g += list_plot(cons[1], plotjoined = True, legend_label = 'Consump2') g += list_plot(Es.flatten(), plotjoined = True, color = 'red', legend_label = 'Control signal') g += list_plot(Y1.flatten(), plotjoined = True, color = 'green', legend_label = 'pressure1') g += list_plot(Y2.flatten(), plotjoined = True, color = 'green', legend_label = 'pressure2') g.show() 
       
##################################################################################################### # This cell implementes the NON-SECURE Pressure control algortihm. # ##################################################################################################### import numpy as np import matplotlib.pylab as plt np.random.seed(1) k = 10 #Period length r = vector(k*[8]) #Pressure requirement T = 5 c = [[0,0,0,0,0,3,3,3,3,3]] #Consumption 1 c1= [[1,1,1,1,1,2,2,2,2,2]] #Consumption 2 cons = np.array([np.repeat(c, T, axis = 0).flatten(), np.repeat(c1, T, axis = 0).flatten()]) E = np.zeros((T, k)) Y1 = np.zeros((T,k)) Y2 = np.zeros((T,k)) kappa = vector(k*[8]) for i in range(T): Y = np.array([kappa - j for j in cons[:,i*k:k*(i+1)]]) Y1[i] = Y[0,:] Y2[i] = Y[1,:] error = np.array([r - j for j in Y]) e_k = vector(np.max(error.transpose(),axis = 1)) E[i] = kappa kappa = kappa + e_k g = Graphics() g += list_plot(cons[0], plotjoined = True, legend_label = 'Consump1') g += list_plot(cons[1], plotjoined = True, legend_label = 'Consump2') g += list_plot(E.flatten(), plotjoined = True, color = 'red', legend_label = 'Control signal') g += list_plot(Y1.flatten(), plotjoined = True, color = 'green', legend_label = 'pressure1') g += list_plot(Y2.flatten(), plotjoined = True, color = 'green', legend_label = 'pressure2') g.show()