adding a communication delay in simulink

I am trying to make a control system block diagram in simulink. I previously made a transfer function with internal communication delay in MATLAB. I am trying to recreate the results in simulink. I used a LTI block to write down in a block, then I simulated the system. Without the delay, I get matching results (with the MATLAB code and simulink) but not with the delay. Is there another way to add the delay?

Answer 1

[Paul]

Try using the LTI System block.

[Kashish Pilyal]

I did use the LTI block ( I mentioned in my question). It did not give the same results as in MATLAB code. I am trying to replicate the results and I am a bit stuck due to this.

[Paul]

I tried the LTI block with a very simple LTI model including a delay and Simuink produced the same output as Matlab. Perhaps there's some mismatch between the Simulink solver settings and the model dynamics. Or maybe something else. Hard to say any more w/o more information.

[Kashish Pilyal]

After the running of simulink model, I linearize it with the model linearizer to get the settling time. The settling time is different.The output graph has a slight overshoot. You mentioned a mismatch between the Simulink solver settings and the model dynamics. How can I check or correct the mismatch?

P.S- I am also using a PD controller in this system

[Paul]

If using a variable step solver, try reducing the Max Step Size parameter.

If using a fixed step solver, try reducing the the Fixed Step Size parameter.

[Kashish Pilyal]

So my transfer function is like where C(s) is the transfer function of PD controller. Can you try a transfer function of this form on your side?

[Paul]

I executed:

s = tf('s');
C = s + 1; theta = 1;
H = 1/(0.5*s + 1)*(exp(-theta*s)*s^2 + C)/(s^2 + C);
step(H)

and then used an LTI Sytem block with H driven by a Step block input. Using the default Simulink solver parameters, the output of the LTI Sytem block overlayed basically perfectly the plot generated by step().

[Kashish Pilyal]

I agree that this does work but if I have the diagram as shown below:

Then this diagram which gives the same tf when made in simulink gives a different answer. I get a difference in the output (the graph has an overshoot). I also get a different settling time. I simulated for delay of 0.02 seconds.

P.S. Sorry the block with should have C(s) as a multiple.

[Paul]

Symbolically, the transfer function from input to output is (assuming I correctly interpeted the P.S. in your comment above)

syms s h theta C(s)
P(s) = 1/(.1*s + 1);
H(s) = (h*s + 1)/s^2;
YoverU(s)= (C(s)/s^2 + exp(-theta*s)) * (P(s)/(1 + P(s)*H(s)*C(s) - 0.8*P(s)))

YoverU(s) = 

[num,den] = numden(YoverU);
YoverU(s) = simplify(num)/simplify(den)

YoverU(s) = 

If you want to get any further here, provide the numerical values for theta and h, and provide C(s) with numbers for the coefficients.

Also, it appears that the LTI System block linearizes delays to a gain of 1. Don't know why it doesn't offer an option to linearize the delays with a Pade approximation.

[Kashish Pilyal]

This image will provide you the data you need. I also simplified and corrected the diagram (sorry for that). I tried using the transport delay block and put 0.02 seconds but it did not work properly. Are there other blocks which I could use?

[Paul]

The new diagram also removed the feedback path with a gain of 0.8.

Anyway, that diagram can be implemented directly in Simulink, except that C(s) will have to be combined with (0.5s + 1)/s^2 into a single Transfer Function block. Is that what you tried?

On the Matlab side, are you trying to get a single LTI system representation from input to output? If so, show your code. If not, please clarify the goal on the Matlab side.

[Kashish Pilyal]

Yes, on Matlab side, I did make a single LTI system (a tf) and gave an input using the lsim command. When I gave the same input to the shown diagram in Simulink, It gave slighlty different results. The output did not match fully as I mentioned there is a slight overshoot in output but if the delay (theta) is set to zero, the results match. This is where I am facing the problem.

[Paul]

I developed the LTI system for that block diagram in Matab and used step()..

In Simulink, I used the LTI system block.

And, in Simulink I impllemented the block diagram explicitly.

All three yielded the same step response.

If you don't show your Matlab code and/or your Simulink diagram, I'm afraid I can't be of any further assistance.

[Kashish Pilyal]

I need to ask that did they all give the same settling times? I do get the same response too but the settling time does not match.

edit: I could send you the models I have if you can confirm this from your side.

[Paul]

All three responses match fine.

You can attach your Simulink model (click the paperclip icon). I won't open it, but maybe someone else will. You can also insert a screen capture of your diagram if you think that would help illustrate what you're doing.

You can copy/paste your Matlab code into a comment here. Preferably formatted as code (clide the code icon to the left of the % icon. If you do that, I can then compare your result to mine based on the block diagram in the figure in this comment

[Kashish Pilyal]

The tf_c and tf_b are:

tf_b=((0.5*s)+1)/(s^2);
tf_c=exp(-0.02*s);

The Matlab code I am using is:

Gamma2=((exp(-0.02*s)*s^2)+(0.6*s)+0.2)/((0.5*s+1)*(s^2+(0.6*s)+0.2));

This should be sufficient for checking if I have done anything wrong. The stepinfo used from the state space from both models gives me a different settling time. Matlab gives 1.82 and the simulink gives 1.95

[Paul]

There are two reasons why the linearized model from Simulink does not match Gamma2.

First, as I already stated in this comment it appears that Simulink does not properly linearize the LTI System block if it includes a delay. Personally, I think this is a bug. So, instead of using the LTI System block for the exp(-theta*s), use a Transport Delay block. In its dialog box, set the Time Delay parameter to theta, and make sure to change the parameter for "Pade order (for linearization)" from the default value of zero. A Pade approximant is basically a transfer function that approximates exp(-s*T). Usually you don't need to set the Pade order higher than 2, but you can play around with it.

Second, the PID Controller block does NOT implement Kp + Ki/s + Kd/s. Rather, it implements Kp + Ki/s + Kd*N/(1 + N/s), where the default value of the filter coefficient paramter, N, is 100. You can see all this in the PID Controller dialog box. So in your derivation of Gamma2, everywhere you had Kd*s, replace it with Kd*100/(1+100/s), assuming you're using the default value of N = 100.

If you really want to use pure PD control of the form of Kp + Kd*s, you can modify the diagram to not use the PID Controller blocks. Delete them both. In the foward path, change the double integrator to a Transfer Fucntion block with Numerator parameter [Kd Kp] and Denominator parameter [1 0 0]. However, by doing this you won't be able to set the integrator initial conditions if you had been doing that. In the feedback path, change tf_b to

tf_b=((0.5*s)+1)/(s^2)*(Kp + Kd*s)

This approach (in addition to using the Transport Delay with Pade approximation) should be nearly identical to Gamma2 after linearization. It won't be exactly identical because of the Pade approximation in the Simulink linearization, but it should be awfully close.

[Kashish Pilyal]

Well thank you for the help but the method you listed gives the same result as the previous one in simulink. The result of Matlab and and simulink settling times is about 0.12 seconds which is close. If this is a bug, is it possible to address it?

[Paul]

I implemented your model in Simulink, using the Transport Delay block with Pade order of 2 and linearized.

>> sys = linearize('pilyal',getlinio('pilyal'));

Then I computed Gamma2, but using the derivative control as it's implemented in the PID block.

>> D = Kd*100/(1+100/s);
>> Gamma2=((exp(-0.02*s)*s^2)+(D)+0.2)/((0.5*s+1)*(s^2+(D)+0.2));

I would have computed it like this:

>> pid = Kp + Kd*100/(1+100/s)
>> clsys2 = series(parallel(exp(-theta*s),pid/s^2),feedback(P,series(pid,H)));

Now compare the stepinfo

>> stepinfo(sys)

ans = 

  struct with fields:

         RiseTime: 1.0609e+00
    TransientTime: 1.8330e+00
     SettlingTime: 1.8330e+00
      SettlingMin: 9.0367e-01
      SettlingMax: 1.0019e+00
        Overshoot: 2.1400e-01
       Undershoot: 0
             Peak: 1.0019e+00
         PeakTime: 3.5920e+00

>> stepinfo(Gamma2)

ans = 

  struct with fields:

         RiseTime: 1.0582e+00
    TransientTime: 1.8185e+00
     SettlingTime: 1.8185e+00
      SettlingMin: 9.0447e-01
      SettlingMax: 1.0022e+00
        Overshoot: 2.2196e-01
       Undershoot: 0
             Peak: 1.0022e+00
         PeakTime: 3.4309e+00
         
>> stepinfo(clsys2)

ans = 

  struct with fields:

         RiseTime: 1.0611e+00
    TransientTime: 1.8339e+00
     SettlingTime: 1.8339e+00
      SettlingMin: 9.0367e-01
      SettlingMax: 1.0019e+00
        Overshoot: 2.0940e-01
       Undershoot: 0
             Peak: 1.0019e+00
         PeakTime: 3.5920e+00
         

Looks pretty close to me.

Also, it turns out there is a way to exactly linearize a model with delays. Search the doc for "Linearize Models with Delays"

[Kashish Pilyal]

I am actually using N=1, for my PID block in simulink. I linearize the model using the model linearizer in simulink. I get the state space and then use the stepinfo. The result is the same as your code of

sys = linearize('pilyal',getlinio('pilyal'));

I checked with both methods. When pade order of 0, the settling time is 1.95 but with pade of 2, I get 1.75 as settling time.

I finally got the settling time of 1.8191 secs. It turns out that the problem was for the PID equation. I was using kp + kd*s while the PID was setting it to . I set the pade approximation to 1 (as I had done the approximation for my thesis before by only taking it until degree 1). When it was set to 0, it did not give proper results (I overlooked the parameter before). Thank you for your help.

In Linearization manager, there is an option to linearize models with delays which gives the settling time of 1.8201, so that is also pretty close.

There are some questions though by this, if it is possible for you to answer them. Why are there slight variations? Is it possible to change the equation in PID block to kp +kd*s instead of ?

[Paul]

When the Pade order is zero, the Tranport Delay block linearizes to a gain of 1. In this case, the closed loop system is simply 1/(0.5*s + 1) regardless of how the PID control C(s) is implemented. in which case we get

stepinfo(tf(1,[.5 1]))

ans = struct with fields:
         RiseTime: 1.0985
    TransientTime: 1.9560
     SettlingTime: 1.9560
      SettlingMin: 0.9045
      SettlingMax: 1.0000
        Overshoot: 0
       Undershoot: 0
             Peak: 1.0000
         PeakTime: 5.2729

The settling time is 1.95, as you've seen.

I changed to N = 1 in my Simulink model with a Transport Delay Pade order of 2 and got a settling time of 1.7979, which is nearly the same result as using the exact delay

s = tf('s');
Kp = 0.1; Kd = 0.6; C = Kp + Kd*s; P = 1/(0.5*s + 1); H = (0.5*s + 1)/s^2; theta = 0.02;
pid1 = Kp + Kd*1/(1+1/s);
% use the exact delay
clsys1 = series(parallel(exp(-theta*s),pid1/s^2),feedback(P,series(pid1,H)));
stepinfo(clsys1)

ans = struct with fields:
         RiseTime: 1.0680
    TransientTime: 1.7972
     SettlingTime: 1.7972
      SettlingMin: 0.9022
      SettlingMax: 1.0048
        Overshoot: 0.4895
       Undershoot: 0
             Peak: 1.0048
         PeakTime: 3.2236

and is the same to four decimal places as replacing the delay with a second order Pade.

% replace the delay with second order Pade
stepinfo(pade(clsys1,2))

ans = struct with fields:
         RiseTime: 1.0682
    TransientTime: 1.7979
     SettlingTime: 1.7979
      SettlingMin: 0.9022
      SettlingMax: 1.0048
        Overshoot: 0.4857
       Undershoot: 0
             Peak: 1.0048
         PeakTime: 3.2236

So with N = 1, I'm seeing a computed settling time of around 1.79, not 1.75.

"Why are there slight variations?"

I'm not sure which variations you're asking about.

"Is it possible to change the equation in PID block to kp +kd*s"

I don't believe so. However, you can construct your own PD using two Gains and a Derivative block. Set the Derivative block parameter "Coefficient c in the transfer function approximation s/(c*s + 1) used for linearization:" to a small number (it has to be positive) to get as close an approximation to pure PD control as you like. This approach might be reasaonable for linearization, but using the Derivative block is not recommended for simulation. Of just use a Transfer Function block fo s/(c*s + 1), keeping in mind that small c means your simulation step size will have be small as well.

Alternatively and perhaps the best approach, as I already metioned in this comment, is to implement the pure PD control by combining it with the double integrator into a single block in the forward path and combining it with tf_b into a single block in the feedback path.

[Kashish Pilyal]

Yes, I did follow the approach as you mentioned. The only difference is if I use the PID block. Is it possible for Matlab devs to do something about it? or I can approach them with this?

[Paul]

If by "do something about it" you mean add an option for the PID block to implement a pure differentiator, you can always submit an enhancement request as described here.

[Kashish Pilyal]

Thank you for the help and I have sent a request.

[Kashish Pilyal]

This is a follow up question to this problem. As I mentioned that I wanted the results of the simulink and the Matlab model to match to some extent. I am facing problem in finding the H infinity norm. For the Matlab model, I get the value of the norm as 1 but for simulink it is infinity. Here is what I am doing.

• I run the simulink model first.
• Then I use the linearization manager to linearize the model at the starting point.
• In the model linearizer box, I use the step function to generate a state space form of the system.
• Then I use the stepinfo command and the hinfnorm command to get the results.

The other parameters match but not the hinfnorm. Is there something wrong that I am doing?

[Paul]

Post the entirety of the code used to develop the "Matab model."

Post a screen capture of the Simulink diagram with any additional information such that I can implement it myself.

[Kashish Pilyal]

The Matlab code is

theta=0.02;
kp=0.2;
kd=0.68;
taup=0.1;
tau=0.1;
h=0.5;
kdd=0;

s=tf('s');
Gamma2=((exp(-theta*s)*s^2)+(kd*s)+kp)/((h*s+1)*(s^2+(kd*s)+kp));

hinfnorm(pade(Gamma2))

The simulink diagram

tf_b=((0.5*s)+1)/(s^2)*(kp + kd*s);

here in the transport delay block, the value of delay is theta and the pade approximant is set to 1.

[Paul]

Repeating the code from above:

theta=0.02;
kp=0.2;
kd=0.68;
taup=0.1;
tau=0.1;
h=0.5;
kdd=0;
s=tf('s');
Gamma2=((exp(-theta*s)*s^2)+(kd*s)+kp)/((h*s+1)*(s^2+(kd*s)+kp));
tf_b=((0.5*s)+1)/(s^2)*(kp + kd*s);

Now sub in the Pade approximant into Gamma2

tempsys = pade(Gamma2)

tempsys =
 
  A = 
          x1     x2     x3     x4     x5     x6     x7     x8     x9
   x1      1      0      0      0      0      0      0      0      0
   x2      0      1      0      0      0      0      0      0      0
   x3      0      0      1      0      0      0      0     -2      0
   x4      0      0      0      1      0      0      0      0      0
   x5      0      0      0      0      1      0      0     -2      0
   x6      0      0      0      0      0  -2.68  -1.56   -0.8      0
   x7      0      0      0      0      0      1      0      0      0
   x8      0      0      0      0      0      0    0.5      0      0
   x9     16      0      0      0      0      0      0      0   -100
 
  B = 
       u1
   x1   0
   x2   0
   x3   0
   x4   0
   x5   0
   x6   2
   x7   0
   x8   0
   x9   0
 
  C = 
         x1    x2    x3    x4    x5    x6    x7    x8    x9
   y1    -1     0     0  0.68     0     0     0   0.4  12.5
 
  D = 
       u1
   y1   0
 
  E = 
       x1  x2  x3  x4  x5  x6  x7  x8  x9
   x1   0   1   0   0   0   0   0   0   0
   x2   0   0   1   0   0   0   0   0   0
   x3   0   0   0   0   0   0   0   0   0
   x4   0   0   0   0   1   0   0   0   0
   x5   0   0   0   0   0   0   0   0   0
   x6   0   0   0   0   0   1   0   0   0
   x7   0   0   0   0   0   0   1   0   0
   x8   0   0   0   0   0   0   0   1   0
   x9   0   0   0   0   0   0   0   0   1
 
Continuous-time state-space model.

We see that tempsys has nine states and is still in descriptor form. Convert to zpk form:

tempsys = zpk(tempsys)

tempsys =
 
  -2 (s-101.4) (s^2 + 0.6748s + 0.1973)
  -------------------------------------
    (s+100) (s+2) (s^2 + 0.68s + 0.2)
 
Continuous-time zero/pole/gain model.

We see that, in reality, tempsys really only has four states in the input/ouput relationship. I'm going to speculate that when presented with a model in descriptor form, hinfnorm() does something similar and gets rid of the extraneous states. The we see that

hinfnorm(tempsys)

ans = 1.0000

On the Simulink side, the linearizer operates in a way that is functionally equivalent to this (I think it actually uses connect() )

D = pade(exp(-theta*s),1);
P = tf(1,[h 1]);
Gamma3 = series(parallel(D,tf([kd kp],[1 0 0 ])),feedback(P,tf_b));

zpk form of Gamma3 shows

zpk(Gamma3)

ans =
 
  -2 s^2 (s-101.4) (s^2 + 0.6748s + 0.1973)
  -----------------------------------------
    s^2 (s+100) (s+2) (s^2 + 0.68s + 0.2)
 
Continuous-time zero/pole/gain model.

Note that Gamma3 has two poles at the origin, which is why

hinfnorm(Gamma3)

ans = Inf

But those poles cancel with the zeros at the origin.

Gamma3 = zpk(minreal(Gamma3))

Gamma3 =
 
  -2 (s-101.4) (s^2 + 0.6748s + 0.1973)
  -------------------------------------
    (s+100) (s+2) (s^2 + 0.68s + 0.2)
 
Continuous-time zero/pole/gain model.

And we see that

hinfnorm(Gamma3)

ans = 1.0000

as expected.

It looks like hinfnorm(), or something it calls, is doing the minimal realization when presented with a descriptor model. Simulink linearizer does not do the minimal realization and leaves that step up to the user (which is the correct thing to do, IMO).

[Kashish Pilyal]

Thank you for the help