**You are now following this question**

- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.

# Pade approximant in Transport delay block

7 views (last 30 days)

Show older comments

##### 0 Comments

### Accepted Answer

Paul
on 22 Jun 2022

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?

##### 25 Comments

Kashish Pilyal
on 22 Jun 2022

Paul
on 22 Jun 2022

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
on 23 Jun 2022

Paul
on 23 Jun 2022

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
on 23 Jun 2022

Paul
on 23 Jun 2022

Edited: Paul
on 23 Jun 2022

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
on 24 Jun 2022

Paul
on 24 Jun 2022

Kashish Pilyal
on 25 Jun 2022

Kashish Pilyal
on 26 Jun 2022

Edited: Kashish Pilyal
on 26 Jun 2022

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
on 26 Jun 2022

Edited: Paul
on 26 Jun 2022

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
on 26 Jun 2022

Paul
on 26 Jun 2022

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
on 26 Jun 2022

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
on 26 Jun 2022

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
on 27 Jun 2022

Edited: Kashish Pilyal
on 27 Jun 2022

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
on 27 Jun 2022

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
on 27 Jun 2022

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
on 27 Jun 2022

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
on 27 Jun 2022

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
on 27 Jun 2022

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
on 27 Jun 2022

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.

Kashish Pilyal
on 28 Jun 2022

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

### More Answers (0)

### See Also

### Categories

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!**An Error Occurred**

Unable to complete the action because of changes made to the page. Reload the page to see its updated state.

Select a Web Site

Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .

You can also select a web site from the following list:

## How to Get Best Site Performance

Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.

### Americas

- América Latina (Español)
- Canada (English)
- United States (English)

### Europe

- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)

- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)

### Asia Pacific

- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)