Convert discrete sys to continuous sys using delays using linear-fractional transformation (LFT)

Hi,

I am trying to create a function that reutrns the continuous time MIMO transfer function of a discrete time sys. Howevee, for some reason although the bodeplot of the continuous system provides the same magnitude magnitude as the discrete sys, the phase angle seems to be mirrored at 0 deg. The function that I would like support is called d2c_exact and is given below.

Could someone figure out what I might be doing wrong? Thanks.

NOTE: I would like a solution where matrices operations are visible and not through matlab functions that obfuscate the process.

% EXAMPLE CODE %

% Declare MIMO matrices A,B,C,D

A = [-0.31261936441776167,0.28772503446055314;0.04760617923328707,0.8887042134058204];

B = [-0.607364659998142,-1.56763974476322,0.5900308503812629;1.4606427684802321,-1.6261530164028424,-0.8222697125090287];

C = [-0.004231925948322231,0.6979760136195045;-0.6958350856345005,-1.5121349753857836];

D = [-1.9523633596883074,-1.4728789255615329,0.14037581452786482;0.22729522014518927,1.2609066299547338,2.830197741550703];

% Define delay of 1 s

Ts = 1;

sys = ss(A,B,C,D,Ts);

% Discrete system

Pd = sys;

% Function that is the topic of the question

Pdc = d2c_exact(sys);

% Matlab equivalent function

Pc = d2c(sys);

% Create plots

wd = logspace(-2,2,200);

opts = bodeoptions;

opts.FreqUnits = 'rad/s';

opts.FreqScale = 'log';

opts.MagUnits = 'abs';

opts.MagScale = 'log';

opts.PhaseUnits = 'deg';

opts.Grid = 'on';

figure

bodeplot(Pd,wd,'b',opts)

hold on

bodeplot(Pdc,wd,'y',opts)

bodeplot(Pc,wd,'r--',opts)

hold off

h = get(gcf,'Children');

xline(h(2),fs/2,'k--');

xline(h(3),fs/2,'k--');

legend('Pd(z)','Pd(s)*inv(ZOH(s))','Pd(z) (ZOH sampling)','Nyquist freq. (fs/2)')

function sysc = d2c_exact(sysd)

nx = size(sysd.A,2);

s=tf('s');

Ts = sysd.Ts;

[A,B,C,D] = ssdata(sysd);

ZOH = exp(-Ts*s); % Transfer function of delay

M = append(ZOH,ZOH);

LR = M - ss(A + eye(nx));

sysc = C*feedback(eye(nx), LR)*B + D;

end

[Paul]

Hi Joan,

The frqeuency response of Pdc is the conjugate of the others.

Do you have a source for the code in d2c_exact?

[Joan]

Hi Paul,

I tried to replicate a code developed in Julia that treats a similar question. The Julia code is in the answer 9/16. Here you can find the link:

Convert discrete to continuous using delays - Specific Domains / Modelling & Simulations - Julia Programming Language (julialang.org)

Answer 1

[Paul]

Hi Joan,

I reviewed the link from the julia site and I think I know what they're doing. But I'm not 100% sure because it seems like they've made the SISO tf solution in post #8 more complicated than it needs to be, and perhaps also in post #9 (which is your link).

As best I can tell, the julia solution can't be implemented in Matlab (I'll show why below), but I think we can do something equivalent.

Here's the code with plots and my comments.

% EXAMPLE CODE %
% Declare MIMO matrices A,B,C,D
A = [-0.31261936441776167,0.28772503446055314;0.04760617923328707,0.8887042134058204];
B = [-0.607364659998142,-1.56763974476322,0.5900308503812629;1.4606427684802321,-1.6261530164028424,-0.8222697125090287];
C = [-0.004231925948322231,0.6979760136195045;-0.6958350856345005,-1.5121349753857836];
D = [-1.9523633596883074,-1.4728789255615329,0.14037581452786482;0.22729522014518927,1.2609066299547338,2.830197741550703];
% Define delay of 1 s
Ts = 1;
sys = ss(A,B,C,D,Ts);
% Discrete system
Pd = sys;
% Function that is the topic of the question
Pdc = d2c_exact(sys);

This isn't really the equivalent function. d2c by default assumes that sys includes the affect of a ZOH on the input, but I don't think that's necessarily an assumption to make here.

% Matlab equivalent function
Pczoh = d2c(sys);

Warning: The model order was increased to handle real negative poles.

For this example, foh works much better

Pcfoh = d2c(sys,'foh');

Warning: The model order was increased to handle real negative poles.

These warning suggests that sys was not developed as discrete-time approxiation using some form of z = exp(s*T).

Here's the new function.

Pdcnew = d2cnew(sys);

Create plots, but only up to the Nyquist frequency, because that's all that makes sense for Pd. Good match, particularly between Pd and Pdcnew.

wd = logspace(-2,2,200);
opts = bodeoptions;
opts.PhaseWrapping = 'on';
bodeplot(Pd,Pdcnew,Pcfoh,wd(wd<pi/Ts),opts)

The model Pdcnew is a system with no states and only internal delays.

Pdcnew

Pdcnew =
 
  D = 
           u1      u2      u3
   y1   7.992   -13.1  -5.396
   y2  -23.19   29.87   15.74
 
  (values computed with all internal delays set to zero)

  Internal delays (seconds): 1  1 
 
Continuous-time state-space model.

Here is your implementation of the post #9 in that julia thread.

function sysc = d2c_exact(sysd)
    nx = size(sysd.A,2);
    s=tf('s');
    Ts = sysd.Ts;
    [A,B,C,D] = ssdata(sysd);

The comment here is correct, i.e., that this is the transfer function of a delay. But the variable name is misleading because this is not the transfer function of a zero-order-hold.

    ZOH = exp(-Ts*s); % Transfer function of delay

But the real problem is that the julia code for that line is

z = delay(-T)

which means that z represents an advance of T seconds (not a delay), which is consistent with the meaning z in a z-transform sense. That means that the equivalent Matlab line would be

z = exp(s*Ts)

But that will cause an error in Matlab because the Control System Toolbox only models delays, not advances.

    M = append(ZOH,ZOH);
    LR = M - ss(A + eye(nx));
    sysc = C*feedback(eye(nx), LR)*B + D;
end

As best I can tell, the julia code approach is to form the block diagram of the discrete-time state space model and then replace z with exp(sT), or equivalently replace z^-1 with exp(-sT), and then form the continuous-time model from the block diagram. I don't understand why they implemented it the way they did, perhaps there's some subtle point that I don't appreciate.

Anyway, here's how that approach (as best I understand it) works in the Control System Toolbox.

function sysc = d2cnew(sysd)

Ts = sysd.Ts;
[A,B,C,D] = ssdata(sysd);

invz = ss(eye(size(A)),'InputDelay',Ts);

%sysc = parallel(series(series(B,feedback(invz,A,+1)),C),D);
sysc = C*feedback(invz,A,+1)*B + D;

end

[Joan]

Hi Paul,

Thanks for you help and support. Your code does the job i.e. the function calculates a continuous equivalent that matches in magnitude and phase in all frequency range (also above the Nyquist frequency). That is what I was trying to achieve.

Thanks for the explanation of the delay. I had the feeling that there was a definition problem there but (as you pointed out) I could not compute the term z = exp(s*Ts) as Matlab's Control Toolbox throw an error. Now I know why. Thanks for that also.

Here is the plot I got with your functions d2cnew:

Answer 2

[Sulaymon Eshkabilov]

Is this what you are getting:

% EXAMPLE CODE %
% Declare MIMO matrices A,B,C,D
A = [-0.31261936441776167,0.28772503446055314;0.04760617923328707,0.8887042134058204];
B = [-0.607364659998142,-1.56763974476322,0.5900308503812629;1.4606427684802321,-1.6261530164028424,-0.8222697125090287];
C = [-0.004231925948322231,0.6979760136195045;-0.6958350856345005,-1.5121349753857836];
D = [-1.9523633596883074,-1.4728789255615329,0.14037581452786482;0.22729522014518927,1.2609066299547338,2.830197741550703];
% Define delay of 1 s
Ts = 1;
sys = ss(A,B,C,D,Ts);
% Discrete system
Pd = sys;
% Function that is the topic of the question
Pdc = d2c_exact(sys);
% Matlab equivalent function
Pc = d2c(sys);

Warning: The model order was increased to handle real negative poles.

% Create plots
wd = logspace(-2,2,200);
% Sampling freq:
fs=1e3;
figure
OPTs = bodeoptions;
OPTs.Grid = 'on';
OPTs.FreqScale = 'linear';
OPTs.Title.String = 'Bode Plot of Transfer Function';
bodeplot(Pd,wd,'b',OPTs)
hold on
bodeplot(Pdc,wd,'y',OPTs)
bodeplot(Pc,wd,'r--',OPTs)
hold off
h = get(gcf,'Children');
xline(h(2),fs/2,'k--');
xline(h(3),fs/2,'k--');
legend('Pd(z)','Pd(s)*inv(ZOH(s))','Pd(z) (ZOH sampling)','Nyquist freq. (fs/2)')

function sysc = d2c_exact(sysd)
    nx = size(sysd.A,2);
    s=tf('s');
    Ts = sysd.Ts;
    [A,B,C,D] = ssdata(sysd);
    ZOH = exp(-Ts*s); % Transfer function of delay
    M = append(ZOH,ZOH);
    LR = M - ss(A + eye(nx));
    sysc = C*feedback(eye(nx), LR)*B + D;
end

[Joan]

Here are the plot options I used. I have edited my answer.

opts = bodeoptions;

opts.FreqUnits = 'rad/s';

opts.FreqScale = 'log';

opts.MagUnits = 'abs';

opts.MagScale = 'log';

opts.PhaseUnits = 'deg';

opts.Grid = 'on';