I want to implement an encrypted controller in Matlab using RSA. I am trying to replicate what is written in this scheme from 'Cyber-Security Enhancement of Networked Control Systems Using Homomorphic Encryption' written by Kiminao Kogiso and Takahiro Fujita:
The input of the controller has to be the output of the plant, and the input of the plant has to be the output of the controller. I think the problem happens when I try to decrypt, but I am not sure.
I am using the state space representation of a system and the following code:
A_p = [9.9998*(10)^-1 1.9700*(10)^-2; -1.970*(10)^-2 9.7025*(10)^-1];
B_p = [9.9900*(10)^-5; 9.8508*(10)^-3];
C_p = [1 0];
D_p = 0;
X0 = [1;0]; %initial conditions for the plant
t = [0:0.01:5]; %sampling time
u = ones(1,length(t)); %encrypted input
y_p = zeros(1, length(t)); %output of the plant
x_p = zeros(2, length(t)); %state of the plant
x_c = zeros(2, length(t)); %state of the controller
x_p(:, 1) = X0; %assign the initial values to the satet of the
%plant
%parameters of the controller
A_c = [1 0.0063; 0 0.3678];
B_c = [0; 0.0063];
C_c = [10 -99.90];
D_c = 3;
end_t = length(t);
N = 94399927; % N=p*q
e = 11; %public key
d = 85800451; %private key
%encrypted parameters for the controller
A_enc = mod(A_c.^e, N);
B_enc = mod(B_c.^e, N);
C_enc = mod(C_c.^e, N);
D_enc = mod(D_c.^e, N);
xc_enc = zeros(2,length(t)); %x_c encrypted
y_enc = zeros(1,length(t)); %y_c encrypted
u_dec = ones(1, length(t)); %u decripted
y_quant = zeros(1, length(t));
for ind=1:end_t-1
%plant
x_p(:, ind+1) = A_p*x_p(: ,ind) + B_p*u_dec(ind);
y_p(ind) = C_p*x_p(: ,ind);
%encryption of x_c and y_p
y_quant(ind) = round(y_p(ind)); %quantization
y_enc(ind) = modpow(y_quant(ind), e, N);
%controller with encrypted parameter
x_c(:, ind+1) = A_enc*x_c(:, ind) - B_enc*y_enc(ind);
u(ind+1) = C_enc*x_c(:, ind) - D_enc*y_enc(ind);
%decrypting u
u_dec(ind+1) = modpow(u(ind+1),d,N);
function result = modpow(base,exp,m)
result = 1;
while (exp > 0)
if bitand(exp,1) > 0
result = mod((result * base),m);
end
exp = bitshift(exp,-1);
base = mod(base^2,m);
end
end
end
figure
plot(t,y_enc);
figure
plot(t,u);
Can somebody help me? Thank's in advance.
EDIT:
sorry, I had to delete some pictures because I can't post all of them because I don't have enough reputation.
Now I am just trying to simply encrypt the output of the plant and then encrypt it before it enters in the controller, and then I encrypt the output of the controller and before it enters in the plant I decrypt it. The idea is the following:
The problem is that if I don't use encryption the outputs I get when I plot
hold on
plot(t, y_p, 'r')
figure
plot(t, u_dec)
are
but when i use encryption, the outputs are
my code is:
A_p = [9.9998*(10)^-1 1.9700*(10)^-2; -1.970*(10)^-2 9.7025*(10)^-1];
B_p = [9.9900*(10)^-5; 9.8508*(10)^-3];
C_p = [1 0];
D_p = 0;
X0 = [1;0];
t = [0:0.01:5];
u = zeros(1,length(t));
y_p = zeros(1, length(t));
x_p = zeros(2, length(t));
x_c = zeros(2, length(t));
x_p(:, 1) = X0;
A_c = [1 0.0063; 0 0.3678];
B_c = [0; 0.0063];
C_c = [10 -99.90];
D_c = 3;
end_t = length(t);
N = 94399927;
e = 11;
d = 85800451;
xc_enc = zeros(2,length(t));
y_enc = zeros(1,length(t));
u_dec = ones(1, length(t));
xc_dec = zeros(2, length(t));
x_dec = zeros(2,length(t));
y_dec = zeros(1, length(t));
u_enc = zeros(1,length(t));
y_quant = zeros(1, length(t));
u_quant = zeros(1, length(t));
for ind=1:end_t-1
x_p(:, ind+1) = A_p*x_p(: ,ind) + B_p*u_dec(ind);
y_p(ind) = C_p*x_p(: ,ind);
y_quant(ind) = round(y_p(ind));
y_enc(ind) = modpow(y_quant(ind),e,N);
y_dec(ind) = modpow(y_enc(ind),d,N);
x_c(:, ind+1) = A_c*x_c(:, ind) - B_c*y_dec(ind);
u(ind+1) = C_c*x_c(:, ind) - D_c*y_dec(ind);
u_quant(ind+1) = round(u(ind+1)); %quantization of u(ind+1)
u_enc(ind+1) = modpow(u_quant(ind+1),e,N);
u_dec(ind+1) = modpow(u_enc(ind+1),d,N);
disp(u(ind +1) + " " + u_dec(ind+1))
end
hold on
plot(t, y_p, 'r')
figure
plot(t, u_dec)
function result = modpow(base,exp,m)
result = 1;
while (exp > 0)
if bitand(exp,1) > 0
result = mod((result * base),m);
end
exp = bitshift(exp,-1);
base = mod(base^2,m);
end
end
function y = powermod_matrix(a,z,n)
% This function calculates y = a^z mod n
%If a is a matrix, it calculates a(j,k)^z mod for every element in a
[ax,ay]=size(a);
% If a is negative, put it back to between 0 and n-1
a=mod(a,n);
% Take care of any cases where the exponent is negative
if (z<0)
z=-z;
for j=1:ax
for k=1:ay
a(j,k)=invmodn(a(j,k),n);
end
end
end
for j=1:ax
for k=1:ay
x=1;
a1=a(j,k);
z1=z;
while (z1 ~= 0)
while (mod(z1,2) ==0)
z1=(z1/2);
a1=mod((a1*a1), n);
end %end while
z1=z1-1;
x=x*a1;
x=mod(x,n);
end
y(j,k)=x;
end %end for k
end %end for j
end
Shouldn't I get almost the same result as the ones without encryption? I think that something is wrong with the input u and how I initialize and increment it, but I don' t know how to fix this. Can somebody help me? Thank's in advance.
EDIT:
moreover, if i add the following line of code:
disp(y_p(ind) + " " + y_dec(ind))
i can see that y_p and y_dec, which should be the same, are different at each iteration.
EDIT:
the problem is that even if i add a quantization ,and so i get integers, i get the same output.
In addition, if i add the following line of code:
disp(u(ind +1) + " " + u_dec(ind+1))
after the decryption of u_enc, so after u_dec, i get :
so u(ind+1) and u_dec(ind+1) are different numbers. And another fact that leaves me confused is that if i quantize u(ind+1), this result doesn't become an integer.
EDIT:
I think that the problem of this is that u(ind+1) becomes a negative number and too big. Infact, if i for example try to do:
x = -19495615862231.37;
z = round(x);
x_enc = modpow(x,11,94399927);
x_dec = modpow(x_enc,85800451,94399927)
as result, after the decryption i get 31779346, and not the original number. I think modpow calculates the modular exponential for large numbers, but it looks like it has a limit.
xc_enc(:, ind+1)
with two different expressions, one after the other. $\endgroup$xc_enc(:, ind+1)
which seems suspicious. I would suggest you to try to see if without encryption your "controller" really controls something. $\endgroup$u_dec
but only update its value at the end of the loop (so the plant is always simulated withu_dec[ind]=1
).y_quant
is never defined.xc_enc
is never used. $\endgroup$