# Pade approximant in Transport delay block

9 次查看（过去 30 天）

显示 更早的评论

### 采纳的回答

Paul
2022-6-22

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 个评论

Kashish Pilyal
2022-6-22

Paul
2022-6-22

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

Paul
2022-6-23

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

Paul
2022-6-23

编辑：Paul
2022-6-23

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

Paul
2022-6-24

Kashish Pilyal
2022-6-25

Kashish Pilyal
2022-6-26

编辑：Kashish Pilyal
2022-6-26

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

编辑：Paul
2022-6-26

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

Paul
2022-6-26

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

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

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

编辑：Kashish Pilyal
2022-6-27

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

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

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

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

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

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

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

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

### 更多回答（0 个）

### 另请参阅

### 类别

### Community Treasure Hunt

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

Start Hunting!**发生错误**

由于页面发生更改，无法完成操作。请重新加载页面以查看其更新后的状态。

您也可以从以下列表中选择网站：

如何获得最佳网站性能

选择中国网站（中文或英文）以获得最佳网站性能。其他 MathWorks 国家/地区网站并未针对您所在位置的访问进行优化。

美洲

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

欧洲

- 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)

亚太

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