Not getting the same step response from Laplace transform and it's respective difference equation

Question

I am trying to simulate a plant on a microcontroller. The transfer function of the plant is

$$ G_{p} \left( s \right) = \frac{2}{\left( s + 3 \right) \left( s - 1 \right)} \tag{1} \label{1}$$

The step response for this function from Octave is

The value goes to $200$ in $6$ seconds and this is what I am trying to reproduce through the difference equation I show a little later.

The Z transform of the above with $T_{s}$ of $0.001 ~ s$ with zero-order hold is

$$ G_{p} \left( z \right) = \frac{9.993 \cdot 10^{-7} z + 9.987 \cdot 10^{-7}}{z^{-2} - 1.998 z + 0.998} \tag{2} \label{2} $$

The difference equation derived from $G_{p} \left( z \right)$ is

$$ y \left( t \right) = 9.993 \cdot 10^{-7} x \left( t - T_{s} \right) + 9.987 \cdot 10^{-7} x \left( t - 2 T_{s} \right) + 1.998 y \left( t - T_{s} \right) - 0.998 y \left( t - 2 T_{s} \right) \tag{3} \label{3} $$

Here is the C code I wrote to realise the above difference equation

#include<stdio.h>

float xtp0 = 0.0;
float etp0 = 0.0;
float xtp0_minus_Ts = 0.0;
float etp0_minus_Ts = 0.0;
float xtp0_minus_2Ts = 0.0;
float etp0_minus_2Ts = 0.0;

float plant0(float input){
    etp0 = input;
    xtp0 =  (9.993e-7F * etp0_minus_Ts)
                    + (9.987e-7F * etp0_minus_2Ts)
                    + (1.998F * xtp0_minus_Ts)
                    - (0.998F * xtp0_minus_2Ts);
    //Saving the history
    xtp0_minus_2Ts = xtp0_minus_Ts;
    etp0_minus_2Ts = etp0_minus_Ts;
    xtp0_minus_Ts = xtp0;
    etp0_minus_Ts = etp0;
    
    
    return xtp0;
}

int main(){
  float x = 0.0F;   
  int i;    
  for(i = 0 ; i < 6000; i++){
     if(i == 0){
       x = plant0(0.0F);     
     }
      else{
     x = plant0(1.0F); 
      }
  } 
  printf("%f\n",x); 
}

I am trying to run the loop $6000$ times as that would amount to $6$ seconds since the sampling period is $0.001$ seconds, and passing the value of $1$ each time I call the plant0 function (thus passing $1$, $6000$ times to plant0). My understanding is that passing a value of $1$ is equivalent to getting the step response of this function. I am expecting the value to be $200$ as observed in the step response graph. However, I get a value of $5.321684$ from the program. Running the same program on the microcontroller is also giving the same output of $5.321684$.

My intention as I stated previously, is to make the difference equation respond in the same way as the step response seen in the plot. Where am I going wrong here?

Answer 1

The difference equation as written in your question is wrong, but I see that you implemented the correct version, also using delayed versions of the output to compute the current output.

The problem is that you used truncated values for the coefficients. You need to represent the coefficients with high accuracy in order for the system to behave as desired. You can't just round to three decimals.

You could actually compute the correct coefficients during an initialization phase in your code, but for a first test you might want to try these more accurate values (b are the numerator coefficients, and a are the denominator coefficients):

b = [9.993339163960613e-07 9.986679158080491e-07]

a = [1 -1.998004995670081e+00 9.980019986673331e-01]

With these coefficient values and with the correct difference equation, the step response of the discrete-time system is indeed just a sampled version of the continuous-time step response, as shown in the figure below: