Pade approximant in Transport delay block

I have some queries about the pade approximant in the transport delay block. I know about the pade approximation and how is it done but I am a little confused as to how it is applied here. Also I assume, if in the transport delay block I increase the value of the order of pade approximant, I will get better results at the cost of longer simulation times. Am I assuming right?

Answer 1

[Paul]

What exactly is the confusion?

Increasing the order of the Pade approximant should have no effect on simulation run time or simulation results. That parameter is only used in developing the linearized model from the simulation block diagram. Have you seen evidence to the contrary?

[Kashish Pilyal]

I was playing around with the values and I noticed that the higher the pade approximant, the higher the states in the state space. When linearizing and finding H infinity norm values, for lower order I had the norm value of 1 and for the higher states as infinity. Also I added a value of 100 to the order for checking. The system became unstable. I thought increasing the value of the order would give better results and not different ones entirely.

[Paul]

Yes, an nth order Pade approximation will add n states to the linearized model

In theory, increasing the order the Pad approximation yields a more accurate approximation of the phase response of the transport delay. But, keep in mind that the Pade approximation is also adding n right-half-plane zeros to your model, which will also have an impact on the results. A 100th order Pade approximation sounds quite excessive and probably will be numerically unreliable, especially if linmod uses a tf representation of the Pade approximation as returned from pade (I don't know that it does and it would be interesting to find out). For example

T = 0.1;
options = bodeoptions;
options.PhaseWrapping = 'on';
options.MagVisible = 'off';
s = tf('s');
h = exp(-s*T);
figure
bodeplot(h,pade(h,2),pade(h,4),logspace(-1,3,500),options)

For a 0.1 second delay, the 2nd order approximant is good to about 30 rad/sec and the 4th order out to about 70-80 rad/sec.

But we run into problems with 100th order approximation

h100 = pade(pade(h,100));
h100.den{:}(end-8:end)

ans = 1×9

1.0e+306 * 0.0001 0.0118 1.6978 Inf Inf Inf Inf Inf Inf

The denominator can't be computed correctly for this value of T!

The bottom line is that the order of the approximation should be as low as possible while still accurately representing the phase of the delay out to some frequency of interest and not causing other numerical problems in the analysis. That frequency of interest will be determined by location of the delay in the system and the properties of the system and/or the properties of the signals being input to the delay.

[Kashish Pilyal]

Ok, but is there an upper limit to the value of pade approximant that I can take? How will I know what is the maximum value that is acceptable? For example take a communication delay of say 10 seconds. In real life scenarios that value is unacceptable but for lower values of that delay, I get a value of H infinity norm of 1 showing that the system is stable.

[Paul]

Now I'm confused.

is there an upper limit to the value of pade approximant

The Pade approximant doesn't have a value, per se. Did you mean an upper limit to the order of the approximant?

but for lower values of that delay

Now are we talking about the delay itself? Or the order of the Pade approximant of that delay?

Is this question really about an upper limit on how much delay the system can tolerate before going unstable?

[Kashish Pilyal]

The question is about an upper limit to the order of approximant before the system goes unstable. I find this puzzling as the I thought that higher the order, the better is the approximation so I would like to know how much higher can I take the value of the order?

[Paul]

Hmm. I've never thought about this question from a theoretical perspective. I'd like to speculate that, in theory, there should be no upper limit on the order of the Pade approximant where instability is induced, just based on a basic Nyquist criterion argument. But that's just speculation on my part.

However, in practice, there is likely a practical upper bound on the order of the approximant before it becomes impractical to use form a computational perspective. As shown above, for T = 0.1, Matlab can't numerically compute the denominator coefficients for N = 100. But problems may arise even for smaller values of N.

T = 0.1;
s = tf('s');
h = exp(-s*T);
pzplot(pade(h,2),pade(h,4),pade(h,20),pade(h,30),ss(pade(h,30)))
legend('N=2','N=4',"N=20","N=30",'N=30 ss')

This plot is attempting to show that the Pade approximant has a pole/zero pattern and that pattern starts to break down somewhere between N = 20 and N = 30 for T = 0.1. Furthermore, the tf form for N=30 can't be reliably converted to ss form, which also indicates concern. As N increases, the poles (and zeros) are migrating further away from origin, the effect of which needs to be understood in the context of any specific analysis.

I'm afraid I can't offer much more on how to determine the upper bound on the order. As stated above, my approach is to use the lowest order such that the phase of the approximant is accurate enough over the frequency range of interest.

[Kashish Pilyal]

It could give me an idea for an approximate upperbound for a given value of a communication delay. As of now, I could set the pade upto 20 order or check which order is stable for a given time delay.

[Paul]

As I tried to say above, increasing the order of the approximant could cause numerical inaccuracies when analyzing the model even before the analysis (likely incorrectly) indicates system instability. Still unclear why such a high order approximation is needed.

[Kashish Pilyal]

The communication delay for lower orders such as 5 or 7, was about 30 seconds and the system was still string stable (H infinity norm was 1) which should not be the case as in real life. Comparing to the other models, I only got proper delay values when pade order was higher but too high would make the system unstable.

[Paul]

Now my curiosity is piqued. Can you post some code of a simple example that illustrates the issue?

[Kashish Pilyal]

So, I think you might have the original system which was in another discussion. I made changes to that.

The other values are as:

tf_a=1/((tau*s)+1);
tf_b=((h*s)+1)/(s^2)*(kp + kd*s + (ki/s));
cons_a=(1-(tau/h));

and the values are:

ki=1e-05;
kp=0.68;
kd=10;
tau=0.1;

The system becomes unstable if the value of H infinity norm is larger than 1. The values of delay and the parameter "h" are varied to get different values of the norm to check that at which values of "h" will we have desirable value of delay such that the system does not become unstable. For the pade order of 2 or 3, the delay can be as big as 40 secs for "h=1.3" while maintaining stability of the system. When you increase the pade order to 20, you get slightly lesser values of delay for which system will be stable (the values for which the value of H infinity norm is equal to 1 and does not go larger than 1).

[Paul]

Have you asked yourself why the system would be unstable for any value of delay and for any order of its Pade approximant? There is no feedback loop around the delay; the delay should have no effect on the stability of the system.

The system becomes unstable if the value of H infinity norm is larger than 1.

Can you clarify that statement? In general, the H-infinity norm is only defined for stable systems and can surely be larger than 1, For example:

h = tf(1,[1 .5 1]);
pole(h)

ans = 
  -0.2500 + 0.9682i
  -0.2500 - 0.9682i

hinfnorm(h)

ans = 2.0656

[Kashish Pilyal]

Yes, it can be but my system is a platoon of cars and if the H infinity norm is greater than 1 then it means that disturbances get passed downstream. It is not just one system but multiple ones. So, if the communicatin delay from a previous vehicle is high in real life, do you really think the controller in the following vehicle would be able to function properly based on information about 30 seconds ago?

[Paul]

An H-inf norm greater than 1 might be undesirable for your application, but it does not indicate instability which is why I found several previous comments above very confusing.

This comment stated: "For the pade order of 2 or 3, the delay can be as big as 40 secs for "h=1.3" while maintaining stability of the system. When you increase the pade order to 20, you get slightly lesser values of delay for which system will be stable (the values for which the value of H infinity norm is equal to 1 and does not go larger than 1)."

Based on where this discussion is now, I think it should have stated: For the pade order of 2 or 3, the delay can be as big as 40 secs for "h=1.3" and keep the H infinity norm of the system less than 1. When you increase the pade order to 20, you get slightly lesser values of delay for which system will have H infinity norm less than 1.

Is my modification to the quoted statement accurate?

If so, I'm still not sure what the issue actually is. Apparently the system with h = 1.3 can tolerate nearly 40 seconds ("as bigs as 40 secs" or "sliightly lesser values") of delay before the H infinity norm exceeds 1. If you're concerned that that result is not sensible, then perhaps you need to reassess whether or not the H-infinity norm is the proper figure-of-merit for this problem.

[Kashish Pilyal]

Well you got the statement accurate but, the H infinity norm does give the maximum value of the output relative to input. So if it is larger than the input (in our case the preceding vehicle), there will be disturbance as the succesor will have more acceleration and it will collide. So I am sure that it is the proper figure-of-merit for this problem.

The main question is why for the values of h=1.3 and delay of 40 (in some cases 50) at pade order of 2-3, I have the norm value of 1 but at pade order of 20 (or some lower values too which I have checked maybe they go upto 5 but not fully sure), and value of "h=1.3" , I get delay of 18 or 19 seconds above which the system becomes unstable. It is a huge difference. Also at this higher order of pade, the maximum delay I can get is approximately 20 seconds and not above that.

[Paul]

Here is my attempt to recreate the results in Matlab for h = 1.3. I also used Simulink w/linearize and saw the same answers.

s = tf('s');
ki=1e-05;
kp=0.68;
kd=10;
tau=0.1;
h = 1.3;
tf_a=1/((tau*s)+1);
tf_b=((h*s)+1)/(s^2)*(kp + kd*s + (ki/s));
cons_a=(1-(tau/h));
pid2 = tf([kd kp ki],[1 0 0 0]);

% 20 second delay
Td = 20;
sys20 = series(parallel(tau/h*exp(-s*Td),pid2),feedback(tf_a,parallel(tf_b,-cons_a)));

% H-inf norm for 5th and 20th order Pade
[hinfnorm(minreal(pade(sys20,5))) hinfnorm(minreal(pade(sys20,20)))]

4 states removed.
4 states removed.

ans = 1×2

1.0000 1.0000

bodemag(sys20,minreal(pade(sys20,5)),minreal(pade(sys20,20)));

4 states removed.
4 states removed.

% 40 second delay
Td = 40;
sys40 = series(parallel(tau/h*exp(-s*Td),pid2),feedback(tf_a,parallel(tf_b,-cons_a)));

% H-inf norm for 5th and 20th order Pade
[hinfnorm(minreal(pade(sys40,5))) hinfnorm(minreal(pade(sys40,20)))]

4 states removed.
4 states removed.

ans = 1×2

1.0000 1.0000

bodemag(sys40,minreal(pade(sys40,5)),minreal(pade(sys40,20)));

4 states removed.
4 states removed.

So this code shows H-inf norm of 1 with Td = 20 or 40 and with Pade order of either 5 or 20. Do these results not match yours?

[Kashish Pilyal]

I am fairly new to the commands of series, parallel and feedback but based on your code of

sys20 = series(parallel(tau/h*exp(-s*Td),pid2),feedback(tf_a,parallel(tf_b,-cons_a)));

The feedback part gave me a picture like this:

which does not match the system, I am aiming for. There is a negative infront of cons_a which is not the case but the tf_b has a negative. Also, I think you forgot the other constant -8.918. If you remove the constant (-8.918), the system should not be stable (H infinity norm greater than 1) for a time delay of more than 20 seconds. So, when you showed the results that it was stable at 40 secs without the constant (-8.918), then I was confused because it should not be the case.

I made changes and

sys20 = series(parallel(tau/h*exp(-s*Td),pid2),feedback(feedback(tf_a,cons_a),-tf_b));

if you run this (which is without the constant) you will get H infinity norm of infinity.

[Paul]

The feedback command assumes negative feedback by default. Because cons_a has positive feedback I negated it on input to parallel. If you don't want to use parallel, the correct command would be

sys20 = series(parallel(tau/h*exp(-s*Td),pid2),feedback(feedback(tf_a,-cons_a),tf_b));

The block diagram in this comment does not include an input perturbation on the signal from the -8.918 constant, and all of the blocks in the model are linear. Consequently, the lineariziation of that block diagram is independent of the -8.918, which is why I ignored it. When you linearize that model the result should be a SISO system from a1input to ai that is independent of the value of the constant block for the disurbance. Is that not what happens?

[Kashish Pilyal]

I am not sure about what happens with the constant block which was leaving me confused. Before, I made the system in simulink, I tried making the system in state space and put the constant in state space but not as an input. I was trying to recreate the same procedure. I would welcome any suggestions on this part.

The state space code gave the same results as the simulink after I linearized the simulink model so I thought that the linearization takes the constant into account.

sys20 = series(parallel(tau/h*exp(-s*Td),pid2),feedback(feedback(tf_a,-cons_a),tf_b));

This code gives the stability even at a delay of 100 secs which seems unfavourable. I also compared the transfer functions that I got from simulink with your code. The ones I have from simulink are different. I am not getting the same results from simulink model when I increase the pade order in transport delay block to 20.

[Paul]

I'm afraid I can't be of any more assistance until you show actual code and results.

Post the complete code that I can copy and paste w/o modification and generate the "state space" code and the sys20. Show Bode (at least the magnitude) plots of both, and on the same plot include the response from the linearized simulink model. Repost a screen shot of the simulink diagram if different than shown above.

Show the results from the hinfnorm command for all three (state space code, sys20, linearized simulink).

[Kashish Pilyal]

h=0.5;
m=-8.918;
tau=0.1;
ki=1e-05;
kp=0.15;
kd=10;

A=[0           0           0         0         0;
   0           0           1         0         0;
   0           0           0         1         0;
   m          -ki         -kp       -kd        0;
   0           0           0        1/h       -1/h];

B=[0;
   0;
   0;
   0;
   1];

C=[0      0         0      -1/h    1/h];

D=[0];
sys=ss(A,B,C,D);
Gamma=tf(sys);
hinfnorm(pade(Gamma))

The code for the state space and the bode plot

For the simulink model, I got a bit different bode plot

You mentioned to show the response, I am confused which one the step response of both systems or what?

[Paul]

The posted code does not appear to match the block diagram in the Simulink model.

1. where is cons_a?

2. where is the time delay?

3. how did m = -8.918 make it into the state space model in the A-matrix? As far as I can tell, -8.8918 doesn't multiply anything the block diagram in the Simulink model. If you want to treat that as a gain on the disturbance, then it should show up in a second column in the B-matrix.

4. why is tau not used?

The result from the ss model is:

h=0.5;
m=-8.918;
tau=0.1;
ki=1e-05;
kp=0.15;
kd=10;
A=[0           0           0         0         0;
   0           0           1         0         0;
   0           0           0         1         0;
   m          -ki         -kp       -kd        0;
   0           0           0        1/h       -1/h];
B=[0;
   0;
   0;
   0;
   1];
C=[0      0         0      -1/h    1/h];
D=[0];
sys=ss(A,B,C,D);
tf(minreal(sys))

4 states removed.

ans =
 
    2
  -----
  s + 2
 
Continuous-time transfer function.

Is that what you expect? Is that also what you get from linearizing the Simulink model?

As best I can tell, that result would only be obtained if cons_a = 0, the transport delay is 0, and tau = h = 0.5 (which would make cons_a = 0), and none of these conditions should be true for the problem at hand, which started out with

tau = 0.1;
h = 1.3;
delay = 20; % or something like that

It's very difficult to provide any assistance if the rules of the game keep changing.

[Kashish Pilyal]

Well to be honest, the system is the same but the equations of the dynamics are a bit different. The one is state space shows the error dynamics while the simulink model just shows the relation ( or a transfer function from the longitudinal dynamics) between the input and the output. If the system is stable or H infinity norm stable, the dynamics no matter which one you use should give the same results. I have those on paper. That is why you will not find tau or cons_a in the state space.

About the delay, I could not put the delay properly in the state space form, that is why I had to use simulink. The step info values from both of the system matched, so after that I decided to add the delay when I started running into problems.

[Paul]

At this point, I really don't have any more to contribute. Good luck with your project.

[Kashish Pilyal]

Ok, thank you for all the assistance. It is really appreciated.