Control system in Matlab with RSA encryption

Question

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.

Answer 1

:EDIT:

Your numbers must be integers POSITIVE INTEGERS. (link)

This is probably the root of the trouble you're having. Interesting that the paper didn't mention anything about this. I'm not sure how they arrived at the encrypted control parameters but, as I mentioned below, it's not my field so I haven't "digested" the paper. In any event, RSA encryption doesn't work with non-integer or negative integer numbers so you can't pass fractional numbers back and forth and you can't pass negative numbers back and forth.

I've removed the original answer here, but it's in the edit history if anyone is curious.

Consider the following example. I used much smaller numbers to reduce my own search times, and I'm sure there are better ways of finding these numbers but, as I mentioned previously, crypto is not my field. Here's some checks:

p = 7
q = 11
N = p*q
lambda = lcm(p-1,q-1)
e = 17
d = 1;
while true
    if(mod(d*e,lambda)==1)
        break;
    else
        d = d+1;
    end
end

>>d

d = 

    23

So now I've got e = 17, d = 23, and N = 77. Now, check that it works:

>> modpow(29,e,N)

ans =

    50

>> modpow(50,d,N)

ans =

    29

Okay, seems to work fine. Now, try it with a negative number, then see if there's a positive number that gives the same encrypted result:

i=1;
while true
    if(modpow(i,e,N)==modpow(-3,e,N))
        break;
    else
        i=i+1;
    end
end

>>i

i = 

    74

Try it again, with a different negative number:

i=1;
while true
    if(modpow(i,e,N)==modpow(-29,e,N))
        break;
    else
        i=i+1;
    end
end

>>i

i =

    48

So, as you can see, a negative integer can wind up giving you the same encryption result as (a different) positive integer.

Since N is always positive (because p and q are prime and thus positive), Matlab documentation states that the modulus operation will always return a positive result. This again confirms that modpow of some negative number, with N>0, will always result in some positive result which could be duplicated by a positive number.

It also means that modpow(~, d, N) will never return a negative number.

Solution

As I mentioned at the top, I'm surprised this wasn't mentioned in the paper you linked, but my approach would be to bias your numbers by some "sufficiently large" positive number and un-bias them before use. Also, since decimal numbers can't be used, you could "inflate" your numbers by multiplying by some scale factor before rounding, then "deflate" the numbers before use by dividing by the same scale factor. Consider the following edits to your code (the powermod_matrix function wasn't being used, so I removed it for brevity):

scaleFactor = 1000; % <-- New
biasValue = 100; % <-- New

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_quant(ind) = round(scaleFactor*(y_p(ind) + biasValue)); % <-- New
     
     y_enc(ind) = modpow(y_quant(ind),e,N);


     y_dec(ind) = modpow(y_enc(ind),d,N);
     
     y_dec(ind) = y_dec(ind)/scaleFactor - biasValue; % <-- New

     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_quant(ind+1) = round(scaleFactor*(u(ind+1) + biasValue)); % <-- New
     
     u_enc(ind+1) = modpow(u_quant(ind+1),e,N);
     u_dec(ind+1) = modpow(u_enc(ind+1),d,N);
     
     u_dec(ind+1) = u_dec(ind+1)/scaleFactor - biasValue; % <-- New

     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

I commented all the lines (six of them) with % <-- New to highlight the changes; again all I did was shift your values by +100 (biasValue) and then multiply by 1000 (scaleFactor) before rounding, and then doing the inverse after decrypting. Seems like that's got your controller working!