adding a communication delay in simulink

Question

Kashish Pilyal 2022-5-26

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1727635-adding-a-communication-delay-in-simulink

评论： Kashish Pilyal 2022-6-13

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?

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Paul 2022-5-27

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1727635-adding-a-communication-delay-in-simulink#answer_972435

Try using the LTI System block.

28 个评论
显示 26更早的评论隐藏 26更早的评论

Paul 2022-5-29

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.

Paul 2022-5-29

编辑：Paul 2022-5-30

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 2022-5-30

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 2022-5-30

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.

Paul 2022-6-12

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 2022-6-13