Implementing digital controller in the time domain

Question

I have simulated a digital control system in the Z domain using MATLAB and I have got satisfactory results. However, when I converted the plant and the digital controller to difference equations and implemented them in an m-file, I did not get the same results displayed in the Z domain

In the Z domain, my system is stable and the response to a unit step met certain desired specifications, but in time domain, the response goes to infinity:

Below you find the program in the Z domain then in the time domain:

1. Z domain simulation:

clear all
close all

J = 0.01;                           % moment of inertia of the rotor       0.01 kg.m^2
b = 0.1;                            % motor viscous friction constant      0.1 N.m.s
K = 0.01;                           % electromotive force constant         0.01 V/rad/sec
                                    % motor torque constant                0.01 N.m/Amp
R = 1;                              % electric resistance                  1 Ohm
L = 0.5;                            % electric inductance                  0.5 H
s = tf('s');
P_motor = K/((J*s+b)*(L*s+R)+K^2)   % Open loop Transfer function of a Dc motor Speed/Voltage
zpk(P_motor)

%rP_motor=0.1/(0.5*s+1)              % First Order approximation of P_motor
Ts = 0.05;
dP_motor = c2d(P_motor, Ts, 'zoh');

zpk(dP_motor)

sys_cl = feedback(dP_motor,1);
[y,t] = step(sys_cl,1);
stairs(t,y);
xlabel('Time (s)')
ylabel('Velocity (rad/s)')
title('Stairstep Response: Original')

Kp = 100;
Ki = 200;
Kd = 10;

C = Kp + Ki/s + Kd*s;
PI=Kp+Ki/s;
PD=Kp+Kd*s;
dC1 = c2d(C,Ts,'tustin')
dPI=c2d(PI,Ts,'tustin');
dPD=c2d(PD,Ts,'tustin');

sys_cl = feedback(dC1*dP_motor,1);
[x2,t] = step(sys_cl,12);
stairs(t,x2)
xlabel('Time (seconds)')
ylabel('Velocity (rad/s)')
title('Stairstep Response: with PID controller')

rlocus(dC1*dP_motor)
axis([-1.5 1.5 -1 1])
title('Root Locus of Compensated System')

z = tf('z',Ts);
dC = dC1/(z+0.82);
rlocus(dC*dP_motor);
axis([-1.5 1.5 -1 1])
title('Root Locus of Compensated System');


sys_cl = feedback(0.8*dC*dP_motor,1);
[x3,t] = step(sys_cl,8);
stairs(t,x3)
xlabel('Time (seconds)')
ylabel('Velocity (rad/s)')
title('Stairstep Response: with Modified PID controller')

Unit step response according to this simulation:

2. Time domain simulation (using difference equations)

close all;

%Output initial condition
yk=0;yk_1=0;yk_2=0;
%Controller initial condition of 
uk_3=0;  uk_2=0;  uk_1=0;  uk=0;          
%Error initial condition
ek_3=0;  ek_2=0;  ek_1=0;  ek=0;  
%Sampling Time
Ts=0.05;

%Desired output (unit step)
yd=1;

%Stokage of different variables (yk, uk and time)
stokuk=[0];
stokyk=[0];
temps=[0];

%******Z Transform of the plant*****%
%     0.002059 z + 0.001686         %
%     ----------------------        %
%     z^2 - 1.511 z + 0.5488        %
%                                   %
%***********************************%

%***Z Transforme of the Controller**%
%    404 z^2 - 632 z + 244          %
%   -------------------------       %
%   z^3 + 0.82 z^2 - z - 0.82       %
%                                   %
%***********************************%



%Calculation of 
for t=[0:Ts:12]
    t;
    
    ek=yd-yk;
    uk=(404*ek_1)-(632*ek_2)+(244*ek_3)-(0.82*uk_1)+(uk_2)+(0.82*uk_3);
    yk=(1.511*yk_1)-(0.5488*yk_2)+(0.002059*uk_1)+(0.001686*uk_2);
    yk_2=yk_1;
    yk_1=yk;
    uk_3=uk_2;
    uk_2=uk_1;
    uk_1=uk;
    ek_3=ek_2;
    ek_2=ek_1;
    ek_1=ek;
    
    stokyk=[stokyk;yk];
    stokuk=[stokuk;uk];
    temps=[temps;t];
     
end

plot(temps,stokyk);grid on
hold on
plot(temps,stokuk);grid on

Unit step response in the time domain:

Could anyone please explain to me why I get this difference in simulations?

Answer 1

Your error is subtle but has a dramatic effect on the stability of the system. Note that we have

$$ e_k = y_d - y_k $$

So you have to compute the error after computing the system output in your loop. If you do it before (like in your code) this is equivalent to adding a time delay of one sample step to your control loop, which in this case is already enough to destroy closed loop stability.

The correct implementation would look like this (copy/pasting most of the code from your question):

close all;

%Output initial condition
yk=0;yk_1=0;yk_2=0;
%Controller initial condition of 
uk_3=0;  uk_2=0;  uk_1=0;  uk=0;          
%Error initial condition
ek_3=0;  ek_2=0;  ek_1=0;  ek=0;  
%Sampling Time
Ts=0.05;

%Desired output (unit step)
yd=1;

%Stokage of different variables (yk, uk and time)
stokuk=[0]; stokyk=[0]; temps=[0];

%Calculation of 
for t=[0:Ts:12]
    t;

    uk=(404*ek_1)-(632*ek_2)+(244*ek_3)-(0.82*uk_1)+(uk_2)+(0.82*uk_3);
    yk=(1.511*yk_1)-(0.5488*yk_2)+(0.002059*uk_1)+(0.001686*uk_2);
    ek=yd-yk; % Compute the control error here, not before yk!

    yk_2=yk_1;
    yk_1=yk;
    uk_3=uk_2;
    uk_2=uk_1;
    uk_1=uk;
    ek_3=ek_2;
    ek_2=ek_1;
    ek_1=ek;
    
    stokyk=[stokyk;yk];
    stokuk=[stokuk;uk];
    temps=[temps;t];
     
end

subplot(2,1,1); stairs(temps,stokyk);grid on;xlabel('t');ylabel('y');ylim([0,1.2]);
subplot(2,1,2); stairs(temps,stokuk);grid on;xlabel('t');ylabel('u');

With this code you get:

which reproduces the step response from your first script.

---

Just one more side note: Having a closed loop system that is so sensitive to such a small increase in the loop delay will likely cause problems in practice. For the sake of robustness you usually want your controller to still be able to stabilize the system even when the delay increases a bit.

You can also see this if you look at the control input $u$ in the figure above: the control signal is oscillating heavily, which is usually undesired.

So if you plan to apply this digital controller in the real world, I would recommend to redesign your controller, with the goal to increase the delay margin of your system.