Convert discrete sys to continuous sys using delays using linear-fractional transformation (LFT)
3 views (last 30 days)
Show older comments
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
2 Comments
Paul
on 20 Dec 2023
Hi Joan,
The frqeuency response of Pdc is the conjugate of the others.
Do you have a source for the code in d2c_exact?
Accepted Answer
Paul
on 20 Dec 2023
Edited: Paul
on 20 Dec 2023
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);
For this example, foh works much better
Pcfoh = d2c(sys,'foh');
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
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
More Answers (1)
Sulaymon Eshkabilov
on 20 Dec 2023
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);
% 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
See Also
Categories
Find more on Uncertainty Analysis in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!