1
$\begingroup$

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

Step response of Gp(s)

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?

$\endgroup$

1 Answer 1

0
$\begingroup$

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:

enter image description here

$\endgroup$
2
  • 1
    $\begingroup$ His denominator is wrong, too. $z^2 - 1.998 z + 0.998$ has poles at $z = 0.998$ and $z = 1$. In translating his one stable, on unstable pole, he's done some rounding and missed the one at $s = +1$. $\endgroup$
    – TimWescott
    Jul 17, 2023 at 19:40
  • $\begingroup$ @TimWescott: You're right. At first I thought that's just the way the coefficients were written in the difference equation, but then I realized that these badly rounded values are also used in the code. $\endgroup$
    – Matt L.
    Jul 17, 2023 at 20:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.