# Control system in Matlab with RSA encryption

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.

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.

• Hi Flavio, I would suggest editing the question to first explain what you are trying, then the code (with maybe some variable clarification) and then clarify what you observe that "doesn't work properly" – N. Staub Dec 4 '18 at 18:30
• Also, there is no control action in your controller and your first control input is arbitrarily equal to 1. I also not that you fill the variable xc_enc(:, ind+1) with two different expressions, one after the other. – N. Staub Dec 4 '18 at 18:36
• you are still updating twice the value xc_enc(:, ind+1) which seems suspicious. I would suggest you to try to see if without encryption your "controller" really controls something. – N. Staub Dec 5 '18 at 8:55
• without encryption i get the correct output i should get if i don't use RSA. The problems come when i encrypt and i can't understand why. – J.D. Dec 5 '18 at 8:58
• You simulated your plant with u_dec but only update its value at the end of the loop (so the plant is always simulated with u_dec[ind]=1). y_quant is never defined. xc_enc  is never used. – N. Staub Dec 5 '18 at 9:11

# :EDIT:

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!

• Thank's for your answer, i have edited the post with the name of the paper i am trying to replicate and at the beginning of the post there is a link to it. – J.D. Dec 6 '18 at 8:03